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In this work we analyze how effects of finite size may modify the thermodynamics of a system of 
strongly interacting fermions that we model using an effective field theory with four-point interac- 
tions at finite temperature and density and look in detail at the case of a confining two-layer system. 
We compute the thermodynamic potential in the large- TV and mean-field approximations and adopt 
a zeta- function regularization scheme to regulate the divergences. Explicit expansions are obtained 
in different regimes of temperature and separation. The analytic structure of the potential is care- 
fully analyzed and relevant integral and series representations for the various expressions involved 
are obtained. Several known results are obtained as limiting case of general results. We numerically 
implement the formalism and compute the thermodynamic potential, the critical temperature and 
the fermion condensate showing that effects of finite size tend to shift the critical points and the 
order of the transitions. The present discussion may be of some relevance for the study of the 
Casimir effect between strongly coupled fermionic materials with inter-layer interactions. 

PACS numbers: 
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I. INTRODUCTION 

The simplest example of an interacting fermion field theory is the Gross- Neveu model in f + f dimensions jl] . Its 
Lagrangian takes the form 

J^ = ^ 7 ^ + ^(#) 2 , (I) 

where N is the number of flavors and g the coupling constant. The model ([1]) enjoys a SU (N) flavor-like symmetry 
and it is invariant under discrete (chiral) Z2-symmetry, ip — > ^ip. Despite its simple form the Gross-Neveu model 
has remarkable properties. It features asymptotic freedom, chiral symmetry breaking and dynamical mass generation, 
making the model a valuable test ground for the more complex case of QCD. Its phase structure, originally studied 
in Ref. [2], is extraordinarily rich and many of its properties were understood only comparatively recently (see, for 
example, Refs. 043)- A very nice review of its phase diagram can be found in Ref. Q- Here, we will limit ourselves 
to recall some of the most salient features of the chiral version of the model, and restrict our considerations to the 
case of homogeneous phases, for which the phase structure can be discussed in terms of the thermodynamic potential 
a i/ . This, using path integrals, can be expressed (we will give details later) as a function of a scalar condensate, 
a ~ (iptp)^ that essentially represents a mass terms for the fermions. Minimizing ^ with respect to a reveals the 
necessary information on the thermodynamics and symmetry breaking patterns of the model. At low temperature, 
for g exceeding a critical value, chiral symmetry is broken and fermions acquire a mass through the appearance of 
a non-vanishing condensate. Once the temperature overcomes a critical value, T c , a phase transition occurs and 
chiral symmetry is restored, signaled by the vanishing of a. Depending on the value of the temperature T and of the 
chemical potential /i, *% has one or two local minima. In the small density - high temperature region, the system is 
in a phase of restored chiral symmetry and has a single minimum at a = 0. When the temperature drops below 
T c , chiral symmetry spontaneously breaks and ^ develops a non trivial minimum. A second order phase transition 
line separates the two phases. The transition changes to first order for larger values of the chemical potential and 
lower values of the temperature. The point where these two lines merge is a tri-critical point, from which two lines 
of metastability depart. Entering (exiting) the region of metastability produces deformations in the thermodynamic 
potential that acquires (loses) a minima. (The above description is valid in the approximation of spatially constant 
condensates. We should mention, however, that inhomogeneous phases are, in fact, favored in some region of the 
phase diagram where the homogeneous approximation becomes too restrictive. Substantial efforts were carried out 
to include such spatially varying phases, and the interested reader may consult Refs. for details and additional 
bibliography. In the following, we will restrict our analysis to the case of spatially constant phases and leave the 
inhomogeneous case for future work. (We will discuss this point in the concluding part of the paper.) 

It is difficult to mention the Gross-Neveu model without referring to its 3 + I-dimensional precursor due to Nambu 
and Jona-Lasinio @ (see Refs. [H, for reviews). In fact, apart from situations in which the relevant spatial 
dimensionality reduces to 1, it seems inevitable to consider extensions of the Gross-Neveu model to 3 + 1 dimensions 
(in 3 + 1-dimensions the model (QJ becomes a special case of the Nambu- Jona Lasinio one with the pseudo-scalar 
contribution suppressed. On the other hand, in 1 + I-dimensions, the model (fT]) augmented by a pseudo-scalar term 
is referred to as NJL2-model, invariant under a continuous chiral symmetry). Extending the Gross-Neveu model to 
higher dimensions is not straightforward due to the fact that such a dimensional lift-up is accompanied by a loss of 
renormalizability. If, at the level of specific applications, this is not a problem per se, and it can be overcome by a 
procedure of phenomenological matching (i.e. fixing the cut-off scale by tuning the model to some experimentally 
measured quantity) , the lack of a natural cut-off scale is often considered as a limitation of fundamental nature. Two 
of the other commonly discussed limitations of the Gross-Neveu model in 1 + 1-dimensions are related to its exact 
integrability, property related to the complete elasticity of the S- matrix [n], [jjj , and to the impossibility of having 
spontaneous symmetry breaking in 1 + 1 dimensions, as stated by the Coleman-Mermin- Wagner theorem [l3], flij . 
Working at lowest order in the approximation of large- N may circumvent the restrictions posed by the Coleman- 
Mermin- Wagner theorem, while dimensional lift-up to 2 + 1-dimensions allows one also to maintain renormalizability 
(see Ref. [15[ for a thorough review) . 

In this paper, we will be concerned with interacting fermion field theories of the form (fT]) in D + 1 dimensions and 
the aspect we aim to analyze is the inclusion of boundaries and finite size effects. Aside of being a natural situation 
to consider in many concrete realistic or semi-realistic applications, it provides a natural possibility to equip these 
models with a physical cut-off scale. 

The simplest example where finite size and boundary effects can be introduced is the n-sphere. Strictly speaking 
this is not a confining cavity, since, in this case, fermions propagate over the surface of the sphere and the effect of 
boundaries is introduced simply via periodicity conditions on the fields. Similar sort of examples can be easily arranged 
by considering toroidal geometries. A number of papers have discussed, in specific contexts, the inclusion of finite size 
effects in the Gross-Neveu and Nambu- Jona Lasinio models on topologically non-trivial geometries and with the fields 
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forced to obey periodic or anti-periodic boundary conditions (see Refs. (16l4l9| ). More precisely, Ref. [16| considered 
the case of the Nambu-Jona Lasinio model at nonzero chemical potential fi and zero temperature on S 1 x S 1 x S 1 
and M. 2 x S 1 and discussed how finite size effects modify the formation of fermion and difermion condensates and the 
generation of a superconducting phase. In Ref. [l7j . the phenomenon of pion condensation was investigated within 
the Nambu-Jona Lasinio model at finite density in M 1 x S 1 . Ref. 18] considered the Gross-Neveu model in three 
dimensions, on § 2 x S 1 and H 2 x S 1 , where § 2 is a 2-sphere and H 2 is a 2-dimensional hyperbolic space. Ref. [l9| 
extended the study of finite size effects within similar class of models to geometries with toroidal and compactified 
directions. See Ref. [20( for earlier results. 

Examples of different sort can be designed by taking a confining enclosure, in which case boundaries arc effectively 
present, and here we will be concerned with the simplest realization of such possibility, namely that of two parallel 
layers. This example, in comparison with those discussed in Refs. fl6l — Tl9l j . presents several complications, as it will 
become clear following the computations reported in the subsequent pages. In fact, while the case of geometries with 
topologically non-trivially circular directions is amenable of a direct analogy with finite temperature quantum field 
theory [HI [22|, for the case of parallel layers the same analogy is much less transparent, due to the non-trivial mixing 
of thermodynamical and geometrical effects. 

Our goal here is to discuss a general method that allows one to include boundary and finite size effects in the 
computation of the thermodynamic potential for a system of interacting fermions. The formalism we will present 
applies to any non-singular geometry, as long as the boundary is smooth, not necessarily connected, and the problem 
separable. A general introduction to the method and a discussion of the assumptions will be presented in Sec. Ill Al 
In Sec. IIIBI we will explicitly perform the computation of the thermodynamic potential for the case of two parallel 
layers with the fermion fields obeying bag boundary conditions at the confining surfaces. Formulae will be given 
for any dimensionality, at non-zero temperature and density, and an explicit expression for the expansion of the 
thermodynamic potential, i.e. the Ginzburg-Landau expansion, for small condensates will be presented up to fourth 
order. The cases of exactly vanishing or large condensates can be treated separately without additional effort. We will 
work in the approximation of mean-field and at leading order in the large- TV limit. Regularization will be performed 
using a zeta- function regularization scheme. Several remarks will be made in Sec. Mil and various limiting cases will 
be considered as check on the general results. Details of the numerical implementation and explicit results will be 
presented in Sec. IIVI for the case of two parallel layers in 3 dimensions. We close the paper with our remarks and a 
brief discussion of a concrete physical system where the present analysis may be of some relevance, namely that of 
the fermionic Casimir effect. 

Results are expressed, as one may expect from the generic form of the effective action, in terms of some series of 
integrals over elliptic functions convoluted with rational functions. Relevant representations and necessary technical 
results are collected in appendices |B] and [C] In Appendix [D] we will briefly discuss how different enclosures can be 
considered within essentially the same formal scheme in a straightforward, although computationally non-trivial, way 
and we will discuss in some details the high-temperature and large condensate cases. 



II. STRONGLY INTERACTING FERMIONS IN BOUNDED ENCLOSURES 



A. Basics 



Consider an interacting fermionic system that can be described by the following action 



dtdv 



D 



9 

2N 



(2) 



where dvu is the volume element in D spatial dimensions. Assume that the system is confined within a region 
^# of D dimensional flat space, bounded by some generic, D — 1 dimensional (co-dimension one) hyper-surface, S, 

where boundary conditions, = are imposed. The boundary £ is assumed to be smooth and may or may not 

? 

be connected. (If the boundary is not connected, e.g. £ = Ei U Sa U ■ • • , the boundary conditions may differ at 
each Ei). Assuming that the Dirac equation is separable (i.e. a sufficient degree of symmetry of the system), one 
may decompose its solutions in positive and negative frequency modes, and express the boundary conditions as a 
quantization condition for the momenta k± perpendicular to the boundary, 



$(|fcj_|) = 0. 



Boundary conditions should be imposed consistently with the requirement of self-adjointness, ensuring the solutions 
to the above eigenvalue equations to be real and positive. This is not the most general case, but simple enough to 
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be analyzed analytically. In general, when it is not possible to find an appropriate coordinate system allowing for 
separation, some sophistication of the approach we present here becomes necessary. More general, non-se par able cases 
can be dealt with using spectral methods, and we will be concerned with this situation somewhere else [23j ]. 
Using the path integral formalism, the partition function can be expressed as a functional integral, 

Q = J @[$, ip] exp (i J drdvD&J , (3) 

where Jz? is the Lagrangian associated with the action ([2]). In the above expression the time direction t has been 
Wick-rotated, t — > it with t E S 1 of radius j3 = 1/T, with T being the temperature and anti-periodicity along 
the S 1 imposed on the fermion fields. It is convenient to re-model the theory by making explicit, by means of a 
Hubbard-Stratonovich transformation, the presence of the condensate, a — —g(\pip)/N, 



n = J @[tp, ip, a] exp^z J drdvD^ef f 



where -Sf e // = ipi'ffid^ip ~ J^ a2 ~ o-ipi/j- The path integral has to be evaluated consistently with the boundary 
conditions imposed on the fields. Sufficiently for our purposes, we can proceed analytically by using the background 
field method and expand the condensate around its classical expectation value, a. Using the large- N and mean-field 
approximations, an expansion of the path integral gives, at leading order, 

p -2 

O = - / dv D ^- + InDet (vy^ - a) + 0(1/ N) , (4) 
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with the trace taken over Dirac and coordinate spaces. The above expression is formal (and independent on the 
boundary conditions that have to be specified at the moment of any concrete computation), and it is valid for 
spatially varying condensates. For the subsequent computation, it is convenient to re-express ((4]) by squaring the 
Dirac operator, 

p -2 i 00 

Q = - / dv D — + - InDet (u>„ - ifif - A + a 2 

I 2q 2 ^-^ L 

J J n=-oo 

+0(1/N), (5) 

where finite density and temperature effects are made explicit and the condensate a restricted to be spatially constant. 
In the above expression ui n = (n + ij are the Matsubara frequencies, (x is a chemical potential and A the Laplacian. 
In order to explicitly compute fi, the use of zeta regularization techniques is particularly advantageous for the present 
case, since it does not require any extra effort when boundary and finite size effects are taken into account. In general, 
the thermodynamic potential can be expressed in terms of the analytically continued values at s — of a generalized 
zeta function, £( s )j an d its derivative (see, for instance, Refs. [HHHI): 

*=i + l^«(o )+m un. ( 6) 

In the present case the function ((s) is associated with the squared Dirac operator appearing in ([5]) and I is a 
renormalization scale. The quantity Vol represents the volume of the geometrical enclosure that confines the system. 
We have discussed details of the zeta-function formalism for interacting fermion field theories in curved spacetimes in 
a previous publication, Ref. [26j j , and the reader is addressed to that paper for details and additional bibliography. 

The requirement of self-adjointness constrains the form of the boundary conditions and, in turn, the form of the 
function $. The simplest (allowed) choice for fermions is that of bag boundary conditions, 

(l + ih a l )ip\z = , (7) 

where n is a unitary vector normal to the boundary . A general discussion on the boundary conditions can be found 
in Ref. [27j . For the case at hand, starting from (J5J), it is straightforward to show that 



r 



■ fci I + uj^ + a 2 — [i 1 — 2ifiUJ Tl 



5 



where D = 3 in the case of three spatial dimensiontQ. Integrating over the direction parallel to the boundaries, one 
obtains 

where the first sum is over the solutions of $(fcj_) = and the second sum is over the Matsubara frequencies. In the 
above expression we have Mellin transformed the summand after integration over k\ \ . The last term in the integrand 
can be expressed in terms of the elliptic function , &2{u,v), using the following identity 

n * 

where 

£pA*) = ^1 + 2 £(-1)" cosh(/3/in)e-^^ (10) 
is, in fact, a series representation for i9 2 ,29]. The sum over fc^ in ([5]) is, instead, recognized as the heat-kernel, 

e(i) = ]Te- tfc i . (ii) 

k ± 

The function &(t) encodes the dependence on the boundary conditions and details of the geometry. Putting things 
together we obtain the following expression 

CM = J" dtf- D ^e-^ 2 e(aH)^ 2 t) , (12) 

where, for convenience, we have rescaled t — > a 2 t. The free energy can be obtained by suitable analytical continuation 
of the above expression and its derivative, according to (j6]), once the set-up is precisely specified. In the next sub- 
section, we will consider the case of two parallel layers. 



B. Two-layers Set-up 

In this section we will focus on the case of parallel flat layers that we assume to be located at z = and z = a, 
leaving the field unconstrained in the remaining directions. For the time being, we will specify the boundary conditions 
to be of the form ([7]). We then may proceed by direct computation. The eigenvalue equation may be expressed, by 
decomposing the solutions of the Dirac equation in positive and negative frequency modes and using ([7| . as an implicit 
constraint for the momenta in the z direction, 

(jsin(fc z a) + k z cos(fc z a) = . (13) 

The roots of the above equation are real and positive, consistently with the self-adjointness of the problem. For 
vanishing er, they can be found explicitly, k z = k = ir(n + 1/2) /a with n € N. (A numerical test on the analytic 
approximations for the eigenvalues developed in this and in the next section is reported in Appendix [A| . 

As it should be clear from (112[) . the first necessary step consists in constructing an explicit expression for the function 
O(i). In general, this cannot be done in fully resummed form. However, in several cases approximate expressions can 
be obtained. Partial knowledge of the phase structure can be obtained in the region of high temperatures, where it 
is possible to use the small-i asymptotics to get a suitable expansion for the thermodynamic potential. Analogously, 
the case of large ct, can also be treated using a similar approximation, since the exponential term in the integrand in 
(fT2|) suppresses the large t portion of the integration range in (fT2|) . However, when finite size effects are taken into 
account, the presence of an additional length scale, the inter-layer separation, complicates things, since the relative 
size of a with respect to a must be taken into account. In the following we will construct the analogous of the 
Ginzburg-Landau-type expansion in all relevant regimes and the remaining part of this subsection will be devoted to 



The case studied in the next section refers to the case of two parallel layers that effectively compactify the direction perpendicular to the 
layers and leave the parallel directions unconstrained. This geometry is clearly non-compact. In this case, we proceed in the usual way 
by closing the cavity with two hyperplanes, separated by a distance L and perpendicular to the boundaries, where we impose periodic 
boundary conditions, and let L to infinity at the end (see, for example, [2£lp. 
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develop an appropriate scheme for this. We will present explicit results up to fourth order, but the procedure can be 
extended to higher orders systematically. 

The two regimes of interest arc those of large aa, relevant to take the infinite volume limit, and that of small aa. 
We will begin considering the latter case. The first step of our scheme consists in solving the quantization condition 
for the momenta k± given by Eq. ([T3|) by perturbatively expanding the eigenvalue equation for small a and fixed a. 
This leads to 

/ fj \ 1 fa\ 2 1 „ /a\ 3 1/1 1 ,\ /ct\ 4 5 / 1 4a 2 \ 

where fco &re the eigenvalues corresponding to a — 0. The above expression, ()14[) . can be used in to arrive at an 
appropriate expansion for 0(t). In the regime of small condensate and small a, we obtain 

Q(a 2 t) = @ Q (a 2 t) - 2(aa)<di(a 2 t) + 2 (aa) 2 G 2 (a 2 t) - - (aa) 3 9 3 (a 2 i) + | (aa) 4 <d 4 (a 2 t) + ■■■ , (15) 

where we have defined 

O (a 2 t) = JtT (t), 

<di(a 2 t) = tJf (t), 
Q 2 (a 2 t)=t 2 ,jeo(t) + ^m(t), 

Q 3 (a 2 t) = t 3 J^(t) + ~i 2 J*i(t) - ~tJti(t) + \tXi(% 



27 15 

e 4 (a 2 t) = t 4 jeo(t) + zt*jent) - 2t 2 je 1 (t) + — t 2 jr 2 (i) - itx 2 (t) + -z-tJtr 3 (t), 



and 



jtr a (t) = Y J (ak )- 2a e- ta k o . 

n=l 

Relation (|15[) can then be used to obtain the following expansion for £(s) 



where 



= 


(4n) D / 2 


(^Ao - 2 (aa) 


A 






Ai 




i , 


A 2 


= ^\s) 




A 3 


= ^\s) 




A 4 




l -24 1 \s) 






n 2s-D 

W "r(.) 



3 v~, -4 



(16) 



2 2 w 2 1 

A 4 = .6/ 4 (0) (s)-2^ 2 (1) (5) + 

with 

dtt s-D/2-l+I e -ta a & p ^( a 2 1) JfTj (t) . (17) 

For convenience we split the integral (|17j) as the sum of two terms, separating out the thermodynamic contribution, 
Wj J \s) from a temperature independent part, Yj J \s), 

^ J) (s) = y} J) (s) + W} J) (s) , (18) 

where 

*f ^ = i^rrE / dtt*-^-^e-^ e -^ , (19) 



m s-DI2-X+I e -ta*<?-^-_^^_ (20) 

Notice that both terms depend on the inter-layer separation that in (j2"0|) is entangled with the temperature. In the 
above expressions we have used the following short-hand notation: v = j + 1/2 and ^2 n v = YlnLi HJLq- Relevant 
representations for the above integrals (fTI?)) and (|2"0"|) and their derivatives along with the analytical continuation 
at s = are constructed in appendix for any dimensionality. Composing the results contained in appendix [B] to 
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obtain the thermodynamic potential presents no special complication, aside of being a tedious and long algebraic 
computation. 

A readable expression for the thermodynamic potential may be presented when the various contributions are 
appropriately combined, 



* = + % + <8i + %T ■ 

The term comes from the ^S J \s) contributions with J = (see appendix IB II part (a)): 



(21) 
(22) 



The coefficients Un and Vn can be computed straightforwardly and for D — 2, 3, 4 are reported in Table (jj) and 
Table Notice that the coefficients ui 1 ^ always occur in finite number and for even spatial dimensionality vanish. 
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TABLE I. Coefficients i4 
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TABLE II. Coefficients v 



(D) 



The J ^ terms can be combined together leading to the following expression ((see appendix IB II part (b))) 



* - ^sbt [£ ^ (■»)" + lo s(i )E tf» («)■ 



4a 



(23) 



where the first few coefficients are tabulated in (|III[) and (|IV|) for D = 2, 3, 4. Finally, the temperature-density 
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TABLE III. Coefficients u 
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contribution, can be expressed as (see Appendix IB 2[) 

2~ r _ _ 2 / 1 

= — — — -r 1 woo - 2 (aer) wio + 2 (acr) w 2 o + -wi; 
+1 (47r)- l V 2 



in terms of the functions 
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TABLE IV. Coefficients v, 



where the quantity on is defined as 

«■ = (l)\i + im 2 . (26) 

The results (|22[) . (|23[) and (p4|) can be added up, according to (|2ip . to obtain the final expression for the thermody- 
namic potential. The above method can be systematically extended to higher orders, although computations become 
increasingly more cumbersome. 



C. Expansion for small a /a 



It is easily noticed that in the regime of separation relatively large compared to the value of the condensate, the 
previous expansion is not useful. When the inter-layer separation is larger than the typical value of tr, the expansion 
developed in the previous section does not converge and to study the thermodynamic behavior in the range where 
a J a is small or aa is large is important to develop an alternative expansion. Clearly, this step is necessary also to 
obtain the infinite volume limit. The aim of this section is to modify the procedure to include the case of large a. 
This case can be done essentially in the same way by taking as expansion parameter a /a. After rescaling the proper 
time t — ¥ a~ 2 t in (IT21 . we expand the function 6(t): 



where 



e(t) = e„(t) 
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(27) 



e (t) 


= ^oW, 
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Using the above expressions we obtain the following expansion for £(s) 



where 
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(s) + ^f)( s ), 



(28) 
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with 



r(a) Jo 



(29) 



Here we will present directly the expression for the thermodynamic potential and skip the details of the computation 
of the above integrals that can be done using a procedure similar to that we have used in appendix for the case of 
small a. 

We express the potential in the following way, 



^9 



(30) 



where + correspond to the sum of the terms with J = 0, "%\ to those with J ^ and to the temperature 
dependent part. Explicitly, + S^o can be expressed in terms of the function 



ri(s) 



T{s + I-D /2) 
T(s) 



2 x D/2-I-s 
2 



(31) 



where the summation can be written in term of the function S*(s, p) defined in Eq. (|B7[) . 
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Proceeding as described in Appendix IB 11 we get 
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The terms and <5^o collect the contributions coming from the first and second terms in curly brackets in the 

expression above, respectively. It takes straightforward steps to arrive at the following formulae: 
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2r-(aa) u + 2r~ (aa) 



D-l 



g r 3 ( acr ) 



D-2 



+ 2 -r-(aa) D -*) log (g) + (r+(aa)^+ 1 - 2r+(aa) D + 2r+{aa) 



D-l 



-r+(aa) 



D-2 



rt (aa) 



D-3 



3 J v ' 3 
where we have adopted the short-hand notation 

D + l 

rls + k 



+ ■ 



= — + r+ + O(s) , for s 







(34) 



(35) 



Notice that the coefficients r^T vanish for even number of spatial dimensions (even D), while only the first (D + l)/2 
are non- vanishing for odd number of spatial dimensions (odd D). We have defined a e = exp ( r y e /2) with 7 e being the 
Euler constant. The contribution <5%> comes from second term of (|B11[) and it can be expressed as 



8% 



D-l 

If 4: 2 ) 



(47T) 



where 



-^(aa) 2 ^ £ q 1 - 2 ^ ^^/^K^o+i (Aqaa) - K,_n^ (2qaa) 

8=1 



(36) 



(37) 
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The J 7^ terms can be combined together leading to the following expression 



2— 1 



V a 



~ ) Tl1 + 2 (a ) V 5 a2ni _ T21 _ Tl2 







-T 2 1 — T 3 l + a T12 — -r 22 — -Tl 3 



(38) 



where 17 j is constructed from the analytical continuation to s = of 

a D+2J T(I- D/2 + s) 



u 



T 2J 



T(s) 



a N D/2-I-s 



(39) 



and its derivative, giving 



s->o \ as 



(40) 



A regular integral representation for the above series can be obtained by use of the Abel-Plana summation formula 
(see appendix B for details): 



tu = ^ ~ # D/2 - / (%+ X- log + a D Vjj 

a 2(I+J) , o_o\fl+£>-2.7-2/V2 / 2 \ 

( s +/ 3 + - Z-j8 log (&r) J 



2 _ 2 \(l+-D-2./-2J)/2 

-a a 



where ^ = ^(0) with 



TT 2 (\ 



and 



Q + iz^j ^ D ' 2 - I (iz) ( X _ log (r 2 ^M) + x+) -(«->■ -2) 



dz 



g27TZ _ ^ ' 



where the constants x±> z± and /3o and /?i are defined in (IC9|) . (|C13|) . (|C14[) and (|C15I) . 

The above integrals can be evaluated analytically by dividing the integration range into two intervals, [0, 1) U [1, oo] 
and by expanding the exponential appropriately. Alternatively, use of numerical approximation does not require any 
effort. 

Finally, the temperature-density contribution, can be expressed as 



%,T = 
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D — 
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a D + 1 (4-K) 
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~3 


-a - 
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-Cl- 











2 ( 






wio + 2 




1 ^20 -1 




.a. 




.a- 







rj^oo - 2 

1 2 3 3 

^30 — 2° ro u + 2^ 21 2, Wl2 



27 



15 



where 



n74o — 2a H72i + 3t373i — 3a n7i2 + — £^22 + -^-^13 



a 2I+2J Cj u 



}■ 



(41) 



(42) 



Combining (|3"4l . (|3^|) . (j3"8")) and (|41l) according to (|28p and (JB]) gives the thermodynamic potential. 
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D. Infinite volume limit 



In the present section we will explicitly take the infinite volume limit. As we shall see, aside for the cancellation of 
spurious divergences, the computations that follow will provide a non-trivial check of the general results. 

To explicitly take the infinite volume limit, rather than the above integral representation we can obtain a simpler 
expression by using the binomial theorem. Choosing aa -C A e N, we have 



u 



,D+2J 



T 2J 



-D-11-ls 



T(I-D/2 + s) 

W) 

v 2 ^ ^fD/2-I-s + k- 1\ ( TT \2fe 
EC" 1 ) ( -/o r . 1 ) (-) Ch(2J -2/,: 1/2) 



D/2-I 




an 



2k-2J 



(43) 



Starting from the above expression, simple power counting arguments are sufficient to show that the infinite volume 
limit can be taken with no problem. The last two terms involve a finite summation over j (y = j + 1/2) and an infinite 
summation over k. Both of them do not produce any diverging behavior, as trivial algebra shows. The first term 
is also regular, since the zeta function occurs with an even argument. Both summations over k are also converging 
in the infinite-volume 'regime', since aa is large. However, using the above representation directly in (|38[) and (|39[) 
requires some care when the limit a — > oo is taken since the leading behavior of the above expression for large a goes 
as STij ~ a 2J and not all terms in ([38| vanish in this limit. Precisely, in the second line of (|38|) the terms proportional 
to a 2 T\2 and — 5T13/2 individually produce linearly diverging contributions as a approaches infinity. These terms, 
when combined, cancel out leaving a regular expression: 



T(s) 



C(4;l/2) C(6;l/2) \ 

7T 4 7T 6 / 



n=0 



n=0 



n=0 



0(l/a) 



+ O(l/o) 



(44) 



Although the representation (|4"3")l is sufficient to prove the regularity of the infinite volume limit, for large, but finite, 
a the integral expression given above is handier for numerical evaluations. 

Letting a — > oo, we have already shown (see ((44]) ) that the second line in ^i, Eq. (f38|) . vanishes. The first and 
fourth terms in the first line of (l38|) also trivially go to zero, while the second and fourth asymptote to constants 
that when combined vanish. The term S^o also vanish in the same limit. The contribution ^.t, is more delicate 
to analyze. All terms in the first line in (|4ip asymptote to zero in the infinite volume limit. In the second line of 
the same expression the first and second terms vanish in the a — > oo limit, while the second and third asymptote to 
a-independent quantities. In the third line of (|4"Tj) . the first, second, third and fifth terms go to zero. The forth and 
sixth terms, while individually diverging as ~ a, when combined, vanish in the a — > oo limit (the reader may check 
this by expanding the quantity — 3n7i2 + 15zj7i3/2 for large a and perform the sums over i exactly). All terms that 
asymptote to a-independent quantities in the infinite volume limit, vanish when the zero temperature limit is also 
taken. Therefore, for T 
giving 



0, only the first term coming from the contribution survives in the a — > oo limit, 



lim % 



(i) 



(D + l)(47r)^ 



a D+1 T 



1 - 



D + l 



(45) 



The above result is seen to reproduce the correct infinite volume zero temperature limit (see Ref. [20| formula (31) 
or Ref. (30j). The terms in < ^ t ,T, that do not vanish in the infinite volume limit, reproduce the correct temperature 
dependence as it can be checked, with some work, by comparing with the flat space limit of the results of Refs. (2(| 
(see also (Hi). 
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III. REMARKS 



In the present section, we wish to make several remarks useful for the subsequent numerical evaluation, as well as 
perform some other checks. 

Vanishing Condensates. The limit of vanishing condensates, a = 0, at zero temperature and chemical potential 
provides one non-trivial check of the result. We will consider the case of D = 3 in order to compare with the correct 
zero temperature free energy for massless free fermions confined between two parallel layers at distance a given in 
Ref. [1H 

This result should be reproduced by our computation in both regimes of small and large separation. In the case of 
small separation, the reader may easily check that the tr-independent contribution comes from the n = term of the 
first summation in (|2"2"|) . This (see Table ((XJ) for the numerical coefficient) reproduces For large separation, the 

leading a-independent term may be computed by expanding (|36[) and by performing exact resummation over q. This 
also reproduces the above result (|46l) . Finally, let us notice that the same result can be obtained directly setting a — 
in (|12[) . In the limit of vanishing temperature and chemical potential, a simple computation gives 

«•) - 2 ^ T± w^ © - » * »> - ' c- ^ °s 

r(s-£>/2) /tt\-2s+£> 



(I) (n(2 S -D)(2^- D -l) 



T(s) \a 

where the function 5? is introduced in (|B7[) . The above expression can be inserted in ©, and reproduces the correct 
result. 



Summations. The results obtained in the preceding sections contain series over certain combinations of Bessel 
functions. Here we will examine the convergence of these series and show how to implement their numerical evalua- 
tion. We will illustrate the procedure for the case D = 3, extensions to any dimensionality being identical. Expression 
(|37jl can be written in terms of the quantity 

oo oo 

£ = J2 & q = ( 22_/ ^-2 (4^) - Kx-i {2qaa)) . (47) 

q=l 9=1 

The above sum may be evaluated in a way that is numerically efficient by separating out the small-q part, 

oo oo 

^ = E^+ E (^-J^l) > (48) 
q=l q=q+i 

with integer q ~ (aa) 1 and the underline signifying expansion for large argument. The first term in the above 
expression can be resummed exactly to any desired order. Each term may contain spurious divergences that can be 
handled by multiplying the summand by q u before resumming over q and then by analytically continuing to u = 0. 
The second term constitute the reminder of the summation and it is expected to be small due to the exponential 
decay for large arguments of the Bessel functions. Computation of the reminder can be handled in a numerically 
efficient way by using the Abel-Plana summation formula (see, for instance, [32l|). 

v 2 - M , n 1 M , \ f°° M , \ , f°° ^{a + iy) ~ J?(a - iy) , 

Yj^^) = 2^(°) + J -^^ dz + 1 J o — — e 4-i — y ■ 

The reminder is found to be negligible in all cases analyzed. 

Numerical Evaluation of the r-functions. Apart from the integral representation for the functions given in the 
preceding section, its possible to evaluate the functions efficiently in a fully numerical way. A convenient approach is 
to fix the values of D, I and J first, then expand tjj(s) and its derivative around s — 0, and then sum the resulting 
expressions numerically. Numerical summations can be performed in a number of ways, here we adopt a sampling 
technique based on computing the summations up to a cut-off, sample the remaining terms, and approximate using 
a polynomial interpolation. The reader may easily check that for I + J > 2 this numerical approach works very 
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efficiently. However, something different is necessary to deal with the case I = J = 1. In this case, we adapt a 
technique similar to that used in Ref. [33| to extend the Chowla-Selberg formula for the Epstein zeta function to 
higher dimensionality. The trick is to write the second term in the summation occurring in (|39|) in terms of its Mellin 
transform, discretize the a direction and integrate over t at each point of the a-grid (the exponential dependence on 
the condensate allows one to interchange the summation and integration operations - see Ref. 33] for details about 
this point in general). The resulting expression can be easily used to construct an interpolated form for tij. 

Convergence. An analogous scheme may be used to compute the function &i t j occurring in the temperature-density 
contribution, however, due to the presence of a chemical potential, the analysis of the convergence of the n summation 

is more delicate. Simple steps show that the series converges exponentially for /i < \J + setting a limit of 

validity for the series representation (|25p . A rather nice consequence of the above inequality is that by decreasing the 
inter-layer separation the convergence improves and allows for larger values of [i than in the corresponding infinite 
volume limit. This is rather convenient for any numerical evaluation, the reason being that for values of a and fi 
satisfying the above inequality the sum over n can be approximated by keeping only the first few terms, and any 
numerical computation does not require any effort. For small /i it is possible to expand the summand appropriately 
and resum each term of the expansion by means of analytical continuation techniques. 

Inter-layer Separation. It is interesting to see that the inter-layer separation, a, appears in the expression for 
the thermodynamic potential in a form identical to that of the Matsubara frequencies (|26[) with an effective tem- 
perature T a = l/(2a). The presence of a chemical potential and higher order terms in the condensate makes the 
combination between geometrical and thermodynamical effects non-trivial, making the scrutiny of the finite-size - 
finite temperature analogy not as transparent as in the case of geometries with topologically non-trivial circular 
directions (22J. 

Scaling. Letting [mass] — > K[mass], where [mass] is any mass parameter (a, fj,, g -1 / 2 , /3 _1 , a -1 and i~ x ), it 
takes simple steps to show that the potential rescales as °i/ — s- k d+1, W , allowing one to fix one of the parameters. 

Large condensates and high temperatures. The high temperature approximation for the potential in the regime 
of large temperature can be obtained with less work than we did above, by using the method of Refs. [H, |35[ which 
consists in rescaling the proper time in the expression for the zeta function (p~2|) . t — > t(3 2 / '(2tt) 2 , and then uses 
the small-i (Schwinger-De Witt) expansion for the heat-kernel. A similar trick can be used to analyze the large 
condensate behavior, as it is easy to understand from the general formula of the zeta function. Rather than reporting 
this computation for the two-layer set-up, which is essentially a repetition of the previous calculation, we will consider 
both cases in some detail in Appendix [Dl for more general boundary geometries. 



IV. RESULTS 



The discussion presented in the previous sub-sections, along with the details reported in the appendices, should 
provide sufficient information on how to proceed with the numerical computation. A technical complication comes 
from having to use different representations in different regimes, but, apart from this, the procedure is straightfor- 
ward. We will present the result of such an evaluation for the 3-dimensional case and 2-dimensional boundaries. 
Generalization to different dimensionality can be performed similarly. 

Zero chemical potential. In the first part of this section we will consider the case of vanishing chemical potential. 
We have already explained the basics of our method in the previous section. One thing that deserves further discussion 
is the way we perform the summations over the Matsubara frequencies and over the eigenvalues on. Clearly, both 
sums can be performed numerically using standard numerical techniques. However, a more efficient way to proceed is 
to expand the Bessel functions in the region of interest and then exactly resum over the Matsubara frequencies each 
term of the resulting series. This provides an analytical expression that will need to be summed over the eigenvalues 
Qij. This second sum is then performed numerically. The full numerical method and partially resummed one are found 
to give consistent results. 

In Fig. [T] we illustrate the result for the thermodynamic potential (normalized by subtracting the value of the 
free energy for a — 0) for sample values of the parameters. Notice that when the inter-layer separation is decreased, the 
effective temperature of the system increases and the minima of the potential shifts towards smaller values. Another 
feature that can be noted is that the potential 'flattens' when the value of a increases. The third prominent feature is 
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that for combinations of distance and temperature it is possible to arrange the system such that the effective potential 
has two minima (one vanishing and one non- vanishing) with equal free energy, suggesting that in certain regions of 
the phase diagram, inhomogeneous phases would appear. Whether even at vanishing density, these could become 
energetically favored owing to finite size effect cannot be proved at this stage, but we will return on this point in the 
future. 




FIG. 1. The figure displays the thermodynamic potential Wn (normalized according to ^| s =o — 0) for vanishing chemical 
potential as a function of the condensate a for fixed inter-layer separation and varying temperature. The values of the parameters 
are fixed to g = 10 2 , 1 -1 = 30. The value of the separation a is indicated in the panels, while the temperature takes the following 
values: = 7. (top curve - left panel), /3 _1 = 1.5 (bottom curve - left panel); = 4.5 (top curve - central panel), = .7 
(bottom curve - central panel); /3 _1 = 2.2 (top curve - right panel), /3 _1 = 0.5 (bottom curve - right panel) 

In Fig. O we illustrate the a-dependence of the potential for fixed values of the temperature. The dashed red line 
represents the critical curve at which a phase transition occurs (for fixed temperature) owing to finite size effects. 

In all figures, increasing the inter-layer separation or decreasing the temperature tends to flatten the potential and 
eventually change the order of the transition from first to second. 




1 2- 3 4 5 I 2 3- 4 5 6 1 2 3 -4 5 6 

a a a 

FIG. 2. The figure illustrates the dependence on the inter-layer separation of the thermodynamic potential for fixed temperature. 
The values of j3 are set as indicated in the panels, while the separation a takes the following values in all panels: a = 5 (top 
curve), a = 20 (bottom curve). The red dashed curves represent the critical curves: (left panel) a — 10.4, (central panel) 
a — 13.2; (right panel) a = 15.12. The other parameters are set as in the previous figure. 

Finite density. Switching on the chemical potential produces some complications related to the fact that its 
dependence modifies the properties of convergence of the expression (|25[). However, as we have stressed previously, 
the incorporation of effects of finite size compensates that of a non- vanishing chemical potential allowing for values of 
/i larger than the corresponding ones in the infinite volume case. The various regimes that can be considered may be 
treated more efficiently using different numerical schemes. For instance, assuming to be in a range of the parameters 
for which the convergence of (|25[) is guaranteed, if \xjT <C 1, it is convenient to expand the summands in (1251) for small 
chemical potentials, proceed by exactly resumming over i and then sum over the Matsubara frequencies by means of 
numerical approximation. When the chemical potential is large compared to the temperature, it is more convenient 
to perform both summations numerically. Here, we always follow this second approach and use the first only for 
check. The double numerical sum is performed over a grid in the tr direction. We perform exactly the sum over the 
two indices up to a cut-off value, sample the additional terms and approximate these using a polynomial fitting. The 
values of the cut-off are varied until the desired accuracy is reached. Results are shown in Fig. ([3]) for sample values 
of the chemical potential and of the inter-layer separation. 
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FIG. 3. The figure illustrates the thermodynamical potential for vanishing vs non-vanishing ^t, fixed a and varying j3. The 
inverse temperature is varied between /3 _1 = 0.6 and /3 _1 = 3 for the left panel, /3 _1 = 1 and = 4.5 for the central 
panel, /3 _1 = 3. and /3 _1 = 7.5 for the right panel. Values of fx and a are indicated in the figures, coupling constant and 
renormalization scale are fixed as in the previous figures. 



The critical temperature, T cr i tl can be computed by solving the equation 

= lim d 2 W/da 2 = 



(49) 



by numerical approximation. We first expand to second order in a the thermodynamic potential, and, after computing 
the second derivative of the resulting expression, we take the limit a — > 0. The infinite summations over the Matsubara 
frequencies and over the eigenvalues on involved are computed using the same numerical sampling technique we have 
used before. The sums are evaluated after discretizing the /3-direction and for fixed values of the chemical potential 
and of the inter-layer separation. The result is then fitted using a polynomial interpolating function and the zero of 
/a./i(/3) = is calculated numerically by using the interpolating function. Alternatively, the critical temperature can 
be computed directly by inspecting how the minima of the potential changes when the temperature and chemical 
potential are varied and by finding for which value of the temperature for fixed fj, and a the non-zero minima of 
the potential vanishes. This method is accurate when the transition is first order while it is less accurate when the 
transition is second order. Results are presented in the right panel of Fig. (j4]). 

We close the section by computing the value that the condensate attains at the minima of the potential, <x m i„, vs 
the temperature (for fixed separation) and vs the separation (for fixed temperature) both at zero chemical potential, 
and the minima are computed using a minimization routine based on the Newton's method. The left panel of Fig. (j4]) 
below shows the value of <J m in obtained for a sample choice of parameters. 
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FIG. 4. The left panel illustrates the behavior of a m in vs the temperature for several values of the inter-layer separation (values 
are indicated in the figure). In the right panel we display the critical temperature, T cr ; t vs /j, for several values of the inter-layer 
separation. The filled dots are directly calculated as explained in the text, while the empty dots are computed by fitting the 
results in the region around the 'knee' and by extrapolating the resulting curve. The renormalization scale and the coupling 
constant are set as in the other figures. 



A precise construction of the three dimensional phase diagram T-\x-a is left for future work as it requires a more 
substantial numerical effort (particularly for the precise determination of the critical points and of the order of the 
transitions). However, the results displayed in Figs. (J2HU) are sufficient to illustrate the main features that the presence 
of the layers produces. As expected, reducing the inter-layer separation induces a shift in the thermodynamic potential 
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indicating an increase in the effective temperature. This is also observed in the tendency of the critical temperature 
to decrease when the separation decreases, indicating that for fixed /3 a decrease of the inter-layer separation tends 
to bring the system towards a phase of restored chiral symmetry. An additional interesting feature is noticed looking 
at the j3 dependence of the potential for fixed a and varying /i (see Fig. ([3])), showing that a simultaneous decrease in 
temperature and separation tends to cause a flattening of the potential and eventually a change in the order of the 
transition from a first order one for large a towards a second order one for small values of a, thus producing possible 
shifts in the critical points. 



V. CONCLUSIONS 

It is not possible, in general, to predict how the finite size of a system may affect its thermodynamical behavior. This 
depends on the geometrical and topological properties of the system, on the presence of external fields, on the specific 
type of interactions, etc. In some cases, symmetries may help to relate high and low temperature-regimes, or (modular 
symmetries) certain topologically non-trivial configurations. Spectral geometric techniques can also be used to make 
more generic statements in specific regimes (e.g. at high temperature). However, a case-by-case analysis seems to 
be necessary to make detailed predictions and understand precisely the combined geometrical and thermodynamical 
effects. 

Here, we have analyzed the impact of finite size effects on the thermodynamics of a system of interacting fermions 
that we have modeled using an effective field theory with a four-fermion interaction term in D dimensions, and 
looked in detail at the case of a confining two-layer system in three dimensions. We have developed a method based 
on zeta-function regularization techniques and used the large-TV and mean-field approximations to construct the 
effective thermodynamic potential. We have developed all necessary expansions in different regimes of temperature 
and separation, constructed all relevant integral and series representations for the various expressions involved and 
thoroughly discussed the analytic structure of the potential. While being particularly interested in the case when the 
'tree-level' four-point interaction term and quantum corrections produced by the deformation in the vacuum due to the 
presence of the boundaries both generate a sizeable force between the boundaries, we have also obtained an expression 
valid for large separations that is important when taking the infinite volume limit. We have performed several checks 
of the individual contributions showing non-trivial cancellation of divergences and reproducing various known results 
in several limits (infinite volume limit, zero condensate, finite temperature). The method we have discussed in this 
paper can be efficiently implemented numerically and, after clarifying several points, we have performed a numerical 
computation of the thermodynamic potential, of the critical temperature, and the condensate in the broken phase for 
the D = 3 dimensional two-layer set-up. 

The most important effect of finite size is the tendency to shift the critical points and eventually the order of the 
transitions. We have seen that the critical temperature and the 'flatness' of the potential depend in a non-trivial 
way on the size of the system (the inter-layer separation in our case) and on the other thermodynamical quantities. 
However, in order to completely describe the modifications that finite size effects may produce to the phase diagram 
of a strongly interacting fermionic system, it is important to relax at least two conditions. The first one is the 
approximation of spatially constant condensates, while the second one is related to the boundary conditions. Relaxing 
the first condition generates several new terms in the partition function and solution of the gap equation in this case 
requires proper functional minimization and related numerical treatment. Regarding the boundary conditions, here 
we have considered the simplest possible choice allowed. However, more general boundary conditions may be used 
[27| . How the phase diagram depends on the boundary conditions remains an open (and in our opinion interesting) 
question and it is clearly very important also in the inhomogeneous case. When the boundaries are disconnected, as 
for instance in the two-layer case analyzed here, different boundary conditions may, in fact, be imposed at each layer 
independently (introducing in the effective potential a dependence on some additional parameter, e.g. a chiral angle). 
While these changes are not expected to modify the situation at high temperature (see for instance appendix [D]) some 
modifications may be expected at low temperature. Especially when inhomogeneous phases become energetically 
favored, how the ground state is affected by the boundary conditions is certainly non-trivial, and whether by means 
of an appropriate choice for the boundary conditions it is possible to engineer the shape of the ground state is also 
an amusing related question. 

A physical situation where the present discussion may have some relevance is in the context of the Casimir effect. 
Here, we have basically considered a Casimir-type set-up with the action augmented by a quartic interaction term. 
An interesting system of this sort is the bilayer graphene with interlayer pairing whose properties have been recently 
discussed in Ref. (36|. (The literature on the fermion Casimir effect is vast. The reader may consult Refs. [37l - [40j 
for some recent articles where the Casimir effect for non-interacting fermions has been discussed.) While in the free 
field case, no phase transition/symmetry breaking occurs, adding an interaction term may change things in a more 
dramatic way. Let us choose the coupling constant g and the temperature in such a way that for large separation a 
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FIG. 5. VEV of the fermion condensate and phases of the Casimir force. The figure shows how the VEV of the condensate 
depends on the inter-layer separation in the region around the transition for fixed temperature (set to = 1.5.). The other 
parameters have been set as in the central panel of Fig. (J2J) - 



chiral symmetry is broken (as in Fig. (j2jl). Starting from this configuration, let us gradually decrease the inter-layer 
separation producing a gradual shift in the minima of the condensate. Once the distance reaches a critical value at 
which chiral symmetry is restored a phase transition occurs: the true minimum of the potential a rn i„ jumps to a 
zero value (see Fig.0 where the transition is a first order one). The transition is from a massive phase (with a mass 
proportional to the condensate) to a massless phase in the restored chiral symmetry region. In the massive phase, the 
fermion contribution to the Casimir force is expected to be suppressed (with the degree of suppression depending on 
the precise value of aa m i n ). On the other hand, in the massless phase, the fermion contribution to the Casimir force 
takes the familiar a -4 behavior. The transition between these two phases of the Casing force (mass-suppressed and 
massless) is related to the spontaneous breaking/restoration of chiral symmetry and it is a phase transition. Fig. ([5]) 
summarizes the above arguments. 

It seems interesting to us to pursue these ideas further and we will return to these issues in the near future [23| ■ 
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Appendix A: Numerical test on the analytical approximation for the eigenvalues 

In order to test the approximations we have developed for the eigenvalues in Sec. Ill Bl and Sec. Ill CI we have 
numerically computed the eigenvalues and compared the resulting values with the analytic approximations over a 
wide range of parameters for k± up to 10 5 . The following figures shows some sample plots of this test for small and 
large values of a and a. The intersection between the vertical red lines and the horizontal axis indicate the eigenvalues 
computed by numerical approximation using a root-finding routine based on the Newton method. The blue dots are 
the eigenvalues obtained by analytical approximation up to fourth order in the expansion parameter. 



massless phase 
Fc ~ a -4 



massive phase 
F c ~ a- 4 e- a "" 
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FIG. 6. Numerical check on the eigenvalue approximation. Plots refer to a = 0.1 and a — 0.1. 
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FIG. 7. Numerical check on the eigenvalue approximation. Plots refer to a = 1 and a — 10. 

Appendix B: Integrals 

In this section we will provide some explicit formulas for the integrals used in the body of the paper (some for- 
mulas and definitions given before will be repeated here for the convenience of the reader). The starting point is the 
following expression 

n 2s-D roo 

^ J \s) = W / dtt*- D / 2 -^e- ta 2ff ^A^tmt) , (Bl) 
1 \ s ) Jo 

where 

.22 

jf j{t ) = n-* J J2^2^> ( B2 ) 

oo 2 2 

■^A a2t ) = l + 2^(-l) n cosh(/3/m)e-^ . (B3) 
n=l 

We have used the short- hand definition v = j + 1/2 and = J^Lo - We wm S P^ t ne integral (IB1[) as 

s/W ( a ) = r 7 (J) (s) + #; (J) (a) , (B4) 



and 



where 



^OO = ^7 ^ S - D/2 - 1+I e- ta2s2 ^r- , (B5) 



X 



dtt s - D / 2 - 1+/ e - ta ff - s ^r t -——. (B6) 
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1. Computation of fj- J \s) 

The first integral can be computed straightforwardly leading to 

rW(s) = a^ni-D/2 + S ) y v _ 2J {a 2- 2+n ^ r - I+ n /2 
i n 2j y(s) ^ y ' 

Let us now consider the two cases J = (a) and J ^ (b) in turn, 
(a) For J — we have 

r (0) (s) = a 2s " Dr ^^ 2 + S ' r ' J ^ 1 ~ 2 -- 2 ' " '~ L ' ' 



T(I 


-D/2 + s) 




T(s) 


T(I 


- D/2 + s) 




T(s) 


T(I 


- D/2 + s) 


T(s) 



a 2s-D 1 v J ~ "I * Z. °> T -2s-2/+n ( » » + ^2 



-S-/+D/2 



, , acr 

-y (s + I-D/2, — 

IT 



2---'- 1 ':/ (s + I-D/2, ^pj 



where we have expressed the sums in terms of the following function 

oo 

y( S ,p):=J2(" 2 + P 2 V S ■ (B7) 

n=l 

The cases a = and a ^ present some slight computational differences. For a — (p = 0), S^{s,p) is nothing but 
the Riemann zeta function, J/*(s,0) = Cr(2s): 

r/ 0) ( S ) = a aa--D £(jjLg^±j) ^-3.-21+1? ( 2 2 s+ 2/-D _ ^ (r (2s + 2J _ D) (Bg) 

r(s) 

For any even spatial dimensionality (any integer and even D), the zeta function is regular and the only poles occur in 
the Gamma functions for s = (denominator) and for s = and I = 0, 1, • • • , -j (numerator) compensating each 
other and leaving a well behaved expression. For any odd spatial dimensionality, any integer and odd D, the zeta 
function has a pole occurring for s = when / = (D + l)/2. This pole is compensated by that of the gamma function 
in the denominator, and the expression is perfectly regular. 

When both a and a are small (that is when p is small), we may expand the function y(s,p) as a series of Riemann 
zeta functions, 

s + q-1 

s 



^(s,p)=^2(-l)"l \p 2 \ R {2 S + 2q) . (B9) 



9=0 

Using this representation we get 

1 ( } r( s ) 1] { s + i-D/2-i 

2q 



x ^2s+2I-D+2 q _ ^ Cfl(2fi + 2/ _ ^ + 2g) 



that recovers (IB8|) in the appropriate limit. As before, for D even and integer, the poles of the Gamma functions 
compensate each other leaving a regular expression for s = 0. For D odd and integer, the poles in the summand, at 
s = and for / = (D + l)/2 — q € N, are canceled by the pole of the Gamma function in the denominator, leaving a 
regular expression. 

A representation valid also for larger values of aa can be constructed by using the Chowla-Selberg representation for 
the series, Ref. (4lj . 



y(s, P ) = + ^^r^P 1 - 2 * + ^f> 2 - S fy- 1/2 ^- 1/2 (2*W) ■ (BIO) 

2 2 r(s) r(s) 
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Using the above expression, we find 

VP (s) = a 2s - D -— — \ y^-T (s + I-D/2- 1/2) (—) + 2tt s+i - d / 2 

1 (SI 2 V 7T / 



(?) 



1/2-S-/+D/2 



V" qS +I-D/2-l/2 
9=1 



2 s + I-D/2 + 1,2 Ks+i _ d/2 _ i/2 (4ga - } 



-K s+I - D/2 -i/2 (2ffaa)] } . (Bll) 

One may notice that the first term in square brackets in (IB10[) cancel out in the above expression. For D odd and 
integer, the poles of the gamma function in in square brackets are compensated by the pole of the gamma function in 
the denominator, leaving a regular expression. For D even integer, the expression vanishes for s = 0. 
(b) For J 7^ 0, we have 

yy\ S ) = a 2s-D W-D/2 + S ) w _ 2J y ^ 2J , 2 _ 2 ^ ,_ J+D/a 



For vanishing a we have 



-2s-2/+£)-2J, -2s-2/+D-2J 



V 



T(s) 

where 



r( s ) 

" \ H (2s + 2/ - D + 2 J; 1/2) , (B12) 



T(I-D/2 + S )__ 2s _ 2I+D _ 2J 



C H (u;b) = J2(k + by u (B13) 



k=0 



is the Hurwitz zeta function. Since the only (simple) pole of the Hurwitz zeta function (|B13I) occur at u = 1, for 
any D even integer the Hurwitz zeta function in (|B12I) is regular for any positive integers / and J, while the poles 
coming from the gamma function in the numerator for I = 0, • • ■ , D/2, are compensated by the poles coming from 
the gamma function in the denominator. For D odd integer, the poles occurring in the Hurwitz zeta function in (IB12|) 
are compensated by those from the gamma function in the denominator. As in the previous cases, the final expression 
is regular. 

When aa is non-vanishing and small, use of the binomial theorem leads to 

y T iJ) (s) = a 2s ~ D F(/ ~ D/2 + S) n -2i+n-2 S -2j (B14) 

1 v ; r(s) v ' 



oo 



q=0 

The regularity of the above expression for s = and for any non zero integer p can be shown as in the previous cases. 
It is also possible to give an integral representation for "fj (s) using, for instance, the Abel-Plana summation formula. 
However, we won't present explicit formulas for this case here. 



2. Computation of #7 <J) (s) 

Wjis) is easier to handle and simple steps lead to 

a 2s — D 00 / /j\ I+s — D/2 

W i J) (S) = ~1W (M) c °sM/3/m) x (B15) 

x n-^Y,^ 2 ' 1 '(*V + a 2 a 2 y W+s - D/2) 'K I+s _ D/2 (n^ 2 V 2 + a 2 a 2 ) 1/2 \ . 

The presence of the Bessel functions in the above expression suppresses, for large values of their argument, the 
summations and make any numerical evaluation straightforward. The regularity is immediate to show and the above 
expression (but not its derivative) vanishes in the limit of s = for any D, I and J. 
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Appendix C: Integral representation of the ^jj-functions. 



Carrying out the analytical continuation of the functions occurring in the regime of large inter-layer separation 
requires an approach slightly different from that we have used in the other cases. The starting point here is the 
function (we suppress the indices / and J for notational convenience and write 3? = J^ij) 



T 2J 



T(s) ^' 



2 \ D/2-I-s 



(CI) 



The case of positive J S N can be dealt with by using some recursive relations starting from the standard Chowla- 
Selberg representation (see 42]). Unfortunately, here we need to consider the case of — J € N and this complicates 
things slightly. Using the Abel-Plana formula, we can express r as 



where 



T (0) = a " F(/ ~ D J 2 + S) 



T 2.1 



T(s) 



T (i) 



a 2J T(I~D/2 + s) 



r 2J 



r( s ) 



v / i \ — 2 J 

i + z ^ D ' 2 - I - S (z)dz 



T (2) 



2) _ a 2J T(I - D/2 + s) 



ir 2J T(s) 



-2J 



\D/2-I-s_ 



-2.1 



where we have used the short-hand notation 
and 













-5 


G- 


-)1 



(z = 0) 



It is easy to get 



where we have defined 



r(o) = ) ^ D/21 [x ~ + s (x+ _ x ~ iog ^ )] + o(s2) ' 

r(7 - D/2 + s) 



T(s) 



X- + S X+ + 0(s 2 ) 



(C2) 



(C3) 



(CM) 



(C5) 



(C6) 



(C7) 



(C8) 



(C9) 



Notice that \- = for D odd. 

In general, the contribution can be expressed in terms of hypergeometric functions. However, a representation 
more useful for our purposes can be constructed in terms of incomplete beta functions, 



p a (x,y) = / t x -\l-t)y-Ht , 



and their regularized form, /3^ eg '(x,y), 



P? a (x,y) = a (x,y)/P{x,v) 



(CIO) 



(Cll) 
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Some computations lead to 



r(1) = Q 2(/+./+s) r(/ _ D/2 + s ) ^_ {1+D -2J-2I-2s)/2 x 



2tt 



r( s ) 



D + l D 

I+J + S — , — -I-a + l I X 



xf3 



(reg) 



I+J + S 



D + l D 



I-s + 1 



where the incomplete gamma function is defined as It takes simple steps to show that the quantity 



r(J - D/2 + s) 
T(s) 



f3 [D/2 - 2(1 - 1 + s), I + J + s - 



D + l 



= + z+s + 0(s 2 ), 



is regular for s — >• 0, while the rest of (|C12j) is regular by construction. Expanding for s —> gives 

a 2(I+J) 



2tt 



-a 2 a*) (1+D - 2J - 2I)/2 x 



x ( Z -pQ + s(z + p Q + Z-p l +2-/3 logCT 2 )) + 0(s 2 ) 



where we have defined 



D + l D 
I + J ^—,-^--1+1 



R _ d f 1+D-2J-2/-3. {reg) 



2 ' 2 



r r D+ 1 D r 

I+J+S _,_-/- s + l 



s=0 



Finally we deal with 



T (2) 



,2.7 



r 2J 



[ x -y + s( X +S -X-J)] + 0(s 2 ) 



where 



j = i 



Q + & D!% -\iz) - (z -> -z) 



dz 



1 

-+ lz 



-2J 



* D I 2 - J (iz) log 0»{iz) - (z -> -z) 



(C12) 



(C13) 

(C14) 
(C15) 

(C16) 



Analytical approximations may be obtained by splitting the integration range into [0, 1) U [1, oo] and by appropriately 
expanding the exponential in the region of small and large z. Both integrals ^ and J? can be computed with no 
effort by means of numerical approximation. 



Appendix D: High temperature and large condensate approximations for general smooth boundaries 

In the present section we wish give some more details on the high-temperature and large-condensate expansions. 
In these cases, it is possible to give simpler and more direct derivations as well as considering slightly more general 
geometries, for the cases where explicit knowledge for the heat-kernel expansion is available. Although it is possible 
to include some singular spaces (for instance conifold geometries, see [43[), here, for simplicity, we will limit our 
discussion to the case of co-dimension one, non-singular boundaries. The starting point, i.e. Eqts. ([5])-([6]), remains 
valid, and the zeta function takes the following general form 

cm = E E (k - + A + s > ( D1 ) 

n X 

where A is used to indicate the eigenvalues of the Laplacian with boundary conditions appropriately imposed. As we 
did before, one may easily recast the above expression in factorized form 

C(s) - r dtt^e-^O (t) (t) , (D2) 

W* 1 I s ) Jo 

where Q(t) = J2\ e ~ tX - 
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1. High Temperature 



The first regime we wish to address is the high temperature one, where we expect the system to be in a region of 
unbroken symmetry. In order to study the high-T regime, we may rescale t — > j^zt, in (|D2j) and arrive at 

2s-l ,.00 / o2 , \ / o2j 



C(*) 



2^T{s) \2irJ J \4tt 2 ) 
At high-temperature (/3 <C 1), we may use the asymptotic expansion of the heat-kernel 



0H 
4^ 



°(£i)=m" D/, s>(£" 

fceN/2 v 



(D3) 



(D4) 



where the coefficients 6^ arc determined by the intrinsic and extrinsic geometry of the spatial section of the background 
manifold as well as by the boundary conditions. The high temperature expansion, skipping lengthy computational 
steps, can be cast as follows 



where, for D = 3. 



W = 
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<P2 



(3) 
¥>3 



(3) 



(D) 



n 
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32tt : 
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2g 16tt 2 
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3/4, v 1/2 = 40F/3, 



-1, 



V3/2 



The coefficients Ufc calculable numbers (for D = 3, the first few are: 

2-v/tt, «2 = 0, • • • ), and we have defined = Vol - 

The above expression is valid for any smooth boundary and for any self-adjoint preserving boundary conditions, 
encoded in the integrated heat-kernel coefficients. Notice that the chemical potential, in the high temperature ex- 
pansion, shows up only at third and higher orders and it multiplies volume (at third order) and boundary (at fourth 
order) terms, but not the condensate (the same is true for bosons 22]). Therefore, up to fourth order in the high 
temperature approximation, it does not affect the position of the minima of the potential. 

The explicit form of the heat-kernel coefficients is not known in general. However, for the present analysis we only 
need information regarding the first few coefficients and this can be retrieved from general formulae (see, for example, 
[3 HE]). In fact, the relevant behavior of the potential can be obtained by looking at the sign of its derivative 
with respect to <7, e = Sign^|^. For the case of bag boundary conditions, the first few coefficients are known and 
sufficient to show that e > for any positive a, implying that the potential is a monotonically increasing function of 
the condensate with an absolute minima at a — 0. Therefore at high temperature and without boundaries chirally 
symmetry is always restored at leading order in the high temperature approximation, as one may have expected. 

For more general boundary conditions (for example, for chiral boundary conditions) the first few coefficients have 
been computed in (45j . The relevant part of 6\ is 



Vol {/ dV TrCr2 + J ^ Tr ( C3Cr757m + C5<7 ^ 5 + c 6°")| 



(D6) 
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where the first integral is over the volume, while the second is over the boundary. The dots indicate terms that either 
vanish or do not depend on a and 7 m indicates the gamma matrices normal to the boundary. The coefficients C3, C5, 
cq are universal multipliers and explicit forms for the case in question can be found in Ref. [45| . Restricting to the 
case of bag boundary conditions that are a special case of chiral boundary condition with the chiral angle set to zero, 
results of Ref. |45( show that C5, cq > and c 3 = 0. Thus, the overall sign of B\ does not change, and the previous 
conclusions remain valid. Changing the boundary conditions allowing for a non-vanishing chiral angle modifies the 
coefficients c. However, for the case of fiat boundaries, the sign of 6\ remains positive. 



2. Large Condensates 



If the values of the condensate is large, one may proceed in a similar way to the high temperature case discussed 
above, owing to the presence of the exponential in (|D2j) that suppresses contributions from the large-i region of the 
integration range. Similar steps give 

l(f <r 3/2 / /3 1/2 p /3 3/2 p 2 

\ (^ 4 + ^r §i/2 ^ " kd * ~ 2 ^ /2 ~ a + ^ k/2 \ + h h + ^' 2 h 

+6 4 ^ - log(ea) 2 Q# a 4 ~ a 2 h + hj) } + y 9 ■ ( D? ) 



where for D = 3 



c = 46*o 5Li 5 /2 

d = 2V2§ 1/2 SLi 2 
15 - 

C2 = tTm^o SLi 7/2 + 26»i 5Li 3/2 



2/3 2 
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T 2 

105.- Jr . 3 
32^° 5Ll »> 2 + W 2 

In the above expression we have defined the following combination of the polylogarithm functions 



c 4 = ^4 ^0 SLig/2 + Tff2~°l SLi 5/2 + ^2 <5^1/2 



SLi a = Li a f- e -^+^ + Li a (-e-^ s -^ (D8) 
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In this work we analyze how effects of finite size may modify the thermodynamics of a system of 
strongly interacting fermions that we model using an effective field theory with four-point interac- 
tions at finite temperature and density and look in detail at the case of a confining two-layer system. 
We compute the thermodynamic potential in the large- TV and mean-field approximations and adopt 
a zeta- function regularization scheme to regulate the divergences. Explicit expansions are obtained 
in different regimes of temperature and separation. The analytic structure of the potential is care- 
fully analyzed and relevant integral and series representations for the various expressions involved 
are obtained. Several known results are obtained as limiting case of general results. We numerically 
implement the formalism and compute the thermodynamic potential, the critical temperature and 
the fermion condensate showing that effects of finite size tend to shift the critical points and the 
order of the transitions. The present discussion may be of some relevance for the study of the 
Casimir effect between strongly coupled fermionic materials with inter-layer interactions. 
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I. INTRODUCTION 

The simplest example of an interacting fermion field theory is the Gross- Neveu model in f + f dimensions jl] . Its 
Lagrangian takes the form 

J^ = ^ 7 ^ + ^(#) 2 , (I) 

where N is the number of flavors and g the coupling constant. The model ([1]) enjoys a SU (N) flavor-like symmetry 
and it is invariant under discrete (chiral) Z2-symmetry, ip — > ^ip. Despite its simple form the Gross-Neveu model 
has remarkable properties. It features asymptotic freedom, chiral symmetry breaking and dynamical mass generation, 
making the model a valuable test ground for the more complex case of QCD. Its phase structure, originally studied 
in Ref. [2], is extraordinarily rich and many of its properties were understood only comparatively recently (see, for 
example, Refs. 043)- A very nice review of its phase diagram can be found in Ref. Q- Here, we will limit ourselves 
to recall some of the most salient features of the chiral version of the model, and restrict our considerations to the 
case of homogeneous phases, for which the phase structure can be discussed in terms of the thermodynamic potential 
a i/ . This, using path integrals, can be expressed (we will give details later) as a function of a scalar condensate, 
a ~ (iptp)^ that essentially represents a mass terms for the fermions. Minimizing ^ with respect to a reveals the 
necessary information on the thermodynamics and symmetry breaking patterns of the model. At low temperature, 
for g exceeding a critical value, chiral symmetry is broken and fermions acquire a mass through the appearance of 
a non-vanishing condensate. Once the temperature overcomes a critical value, T c , a phase transition occurs and 
chiral symmetry is restored, signaled by the vanishing of a. Depending on the value of the temperature T and of the 
chemical potential /i, *% has one or two local minima. In the small density - high temperature region, the system is 
in a phase of restored chiral symmetry and has a single minimum at a = 0. When the temperature drops below 
T c , chiral symmetry spontaneously breaks and ^ develops a non trivial minimum. A second order phase transition 
line separates the two phases. The transition changes to first order for larger values of the chemical potential and 
lower values of the temperature. The point where these two lines merge is a tri-critical point, from which two lines 
of metastability depart. Entering (exiting) the region of metastability produces deformations in the thermodynamic 
potential that acquires (loses) a minima. (The above description is valid in the approximation of spatially constant 
condensates. We should mention, however, that inhomogeneous phases are, in fact, favored in some region of the 
phase diagram where the homogeneous approximation becomes too restrictive. Substantial efforts were carried out 
to include such spatially varying phases, and the interested reader may consult Refs. for details and additional 
bibliography. In the following, we will restrict our analysis to the case of spatially constant phases and leave the 
inhomogeneous case for future work. (We will discuss this point in the concluding part of the paper.) 

It is difficult to mention the Gross-Neveu model without referring to its 3 + 1-dimensional precursor, the Nambu- 
Jona Lasinio (NJL) model [8| (see Refs. for reviews). In fact, apart from situations in which the relevant spatial 

dimensionality reduces to 1, it seems inevitable to consider extensions of the Gross-Neveu model to 3 + 1 dimensions 
(in 3 + I-dimensions the model (QJ becomes a special case of the Nambu-Jona Lasinio one with the pseudo-scalar 
contribution suppressed. On the other hand, in I + I-dimensions, the model (fT]) augmented by a pseudo-scalar term 
is referred to as NJL2-model, invariant under a continuous chiral symmetry). Extending the Gross-Neveu model to 
higher dimensions is not straightforward due to the fact that such a dimensional lift-up is accompanied by a loss of 
renormalizability. If, at the level of specific applications, this is not a problem per se, and it can be overcome by a 
procedure of phenomenological matching (i.e. fixing the cut-off scale by tuning the model to some experimentally 
measured quantity) , the lack of a natural cut-off scale is often considered as a limitation of fundamental nature. Two 
of the other commonly discussed limitations of the Gross-Neveu model in I + 1-dimensions are related to its exact 
integrability, property related to the complete elasticity of the S- matrix [n], [jjj , and to the impossibility of having 
spontaneous symmetry breaking in I + I dimensions, as stated by the Coleman-Mermin- Wagner theorem [l3], flij . 
Working at lowest order in the approximation of large- N may circumvent the restrictions posed by the Coleman- 
Mermin- Wagner theorem, while dimensional lift-up to 2 + I-dimensions allows one also to maintain renormalizability 
(see Ref. [15[ for a thorough review) . 

In this paper, we will be concerned with interacting fermion field theories of the form ([T]) in D + 1 dimensions and 
the aspect we aim to analyze is the inclusion of boundaries and finite size effects. Aside of being a natural situation 
to consider in many concrete realistic or semi-realistic applications, it provides a natural possibility to equip these 
models with a physical cut-off scale. 

The simplest example where finite size and boundary effects can be introduced is the n-sphere. Strictly speaking 
this is not a confining cavity, since, in this case, fermions propagate over the surface of the sphere and the effect of 
boundaries is introduced simply via periodicity conditions on the fields. Similar sort of examples can be easily arranged 
by considering toroidal geometries. A number of papers have discussed, in specific contexts, the inclusion of finite size 
effects in the Gross-Neveu and Nambu-Jona Lasinio models on topologically non-trivial geometries and with the fields 
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forced to obey periodic or anti-periodic boundary conditions (see Refs. |16l4l9l |). More precisely, Ref. [16| considered 
the case of the Nambu-Jona Lasinio model at nonzero chemical potential fi and zero temperature on S 1 x S 1 x S 1 
and M. 2 x S 1 and discussed how finite size effects modify the formation of fermion and difermion condensates and the 
generation of a superconducting phase. In Ref. [l7j . the phenomenon of pion condensation was investigated within 
the Nambu-Jona Lasinio model at finite density in M 1 x S 1 . Ref. 18] considered the Gross-Neveu model in three 
dimensions, on § 2 x S 1 and H 2 x S 1 , where § 2 is a 2-sphere and H 2 is a 2-dimensional hyperbolic space. Ref. [l9| 
extended the study of finite size effects within similar class of models to geometries with toroidal and compactified 
directions. See Ref. [20( for earlier results. 

Examples of different sort can be designed by taking a confining enclosure, in which case boundaries arc effectively 
present, and here we will be concerned with the simplest realization of such possibility, namely that of two parallel 
layers. This example, in comparison with those discussed in Refs. fl6l — Tl9l j . presents several complications, as it will 
become clear following the computations reported in the subsequent pages. In fact, while the case of geometries with 
topologically non-trivially circular directions is amenable of a direct analogy with finite temperature quantum field 
theory [HI [22|, for the case of parallel layers the same analogy is much less transparent, due to the non-trivial mixing 
of thermodynamical and geometrical effects. 

Our goal here is to discuss a general method that allows one to include boundary and finite size effects in the 
computation of the thermodynamic potential for a system of interacting fermions. The formalism we will present 
applies to any non-singular geometry, as long as the boundary is smooth, not necessarily connected, and the problem 
separable. A general introduction to the method and a discussion of the assumptions will be presented in Sec. Ill Al 
In Sec. IIIBI we will explicitly perform the computation of the thermodynamic potential for the case of two parallel 
layers with the fermion fields obeying bag boundary conditions at the confining surfaces. Formulae will be given 
for any dimensionality, at non-zero temperature and density, and an explicit expression for the expansion of the 
thermodynamic potential, i.e. the Ginzburg-Landau expansion, for small condensates will be presented up to fourth 
order. The cases of exactly vanishing or large condensates can be treated separately without additional effort. We will 
work in the approximation of mean-field and at leading order in the large- TV limit. Regularization will be performed 
using a zeta- function regularization scheme. Several remarks will be made in Sec. Mil and various limiting cases will 
be considered as check on the general results. Details of the numerical implementation and explicit results will be 
presented in Sec. IIVI for the case of two parallel layers in 3 dimensions. We close the paper with our remarks and a 
brief discussion of a concrete physical system where the present analysis may be of some relevance, namely that of 
the fermionic Casimir effect. 

Results are expressed, as one may expect from the generic form of the effective action, in terms of some series of 
integrals over elliptic functions convoluted with rational functions. Relevant representations and necessary technical 
results are collected in appendices |B] and |DJ In Appendix [E] we will briefly discuss how different enclosures can be 
considered within essentially the same formal scheme in a straightforward, although computationally non-trivial, way 
and we will discuss in some details the high-temperature and large condensate cases. 



II. STRONGLY INTERACTING FERMIONS IN BOUNDED ENCLOSURES 



A. Basics 



Consider an interacting fermionic system that can be described by the following action 



dtdv 



D 



9 

2N 



(2) 



where dvu is the volume element in D spatial dimensions. Assume that the system is confined within a region 
^# of D dimensional flat space, bounded by some generic, D — 1 dimensional (co-dimension one) hyper-surface, S, 

where boundary conditions, = are imposed. The boundary £ is assumed to be smooth and may or may not 

? 

be connected. (If the boundary is not connected, e.g. £ = Ei U Sa U ■ • • , the boundary conditions may differ at 
each Ei). Assuming that the Dirac equation is separable (i.e. a sufficient degree of symmetry of the system), one 
may decompose its solutions in positive and negative frequency modes, and express the boundary conditions as a 
quantization condition for the momenta k± perpendicular to the boundary, 



$(|fcj_|) = 0. 



Boundary conditions should be imposed consistently with the requirement of self-adjointness, ensuring the solutions 
to the above eigenvalue equations to be real and positive. This is not the most general case, but simple enough to 



4 



be analyzed analytically. In general, when it is not possible to find an appropriate coordinate system allowing for 
separation, some sophistication of the approach we present here becomes necessary. More general, non-se par able cases 
can be dealt with using spectral methods, and we will be concerned with this situation somewhere else [23j ]. 
Using the path integral formalism, the partition function can be expressed as a functional integral, 

Q = J @[$, ip] exp (i J drdvD&J , (3) 

where Jz? is the Lagrangian associated with the action ([2]). In the above expression the time direction t has been 
Wick-rotated, t — > it with t E S 1 of radius j3 = 1/T, with T being the temperature and anti-periodicity along 
the S 1 imposed on the fermion fields. It is convenient to re-model the theory by making explicit, by means of a 
Hubbard-Stratonovich transformation, the presence of the condensate, a — —g(\pip)/N, 



n = J @[tp, ip, a] exp^z J drdvD^ef f 



where -Sf e // = ipi'ffid^ip ~ J^ a2 ~ o-ipi/j- The path integral has to be evaluated consistently with the boundary 
conditions imposed on the fields. Sufficiently for our purposes, we can proceed analytically by using the background 
field method and expand the condensate around its classical expectation value, a. Using the large- N and mean-field 
approximations, an expansion of the path integral gives, at leading order, 

p -2 

O = - / dv D ^- + InDet (vy^ - a) + 0(1/ N) , (4) 
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with the trace taken over Dirac and coordinate spaces. The above expression is formal (and independent on the 
boundary conditions that have to be specified at the moment of any concrete computation), and it is valid for 
spatially varying condensates. For the subsequent computation, it is convenient to re-express ((4]) by squaring the 
Dirac operator, 

p -2 i 00 

Q = - / dv D — + - InDet (u>„ - ifif - A + a 2 

I 2q 2 ^-^ L 

J J n=-oo 

+0(1/N), (5) 

where finite density and temperature effects are made explicit and the condensate a restricted to be spatially constant. 
In the above expression ui n = (n + ij are the Matsubara frequencies, (x is a chemical potential and A the Laplacian. 
In order to explicitly compute fi, the use of zeta regularization techniques is particularly advantageous for the present 
case, since it does not require any extra effort when boundary and finite size effects are taken into account. In general, 
the thermodynamic potential can be expressed in terms of the analytically continued values at s — of a generalized 
zeta function, £( s )j an d its derivative (see, for instance, Refs. [HHHI): 

*=i + l^«(o )+m un. ( 6) 

In the present case the function ((s) is associated with the squared Dirac operator appearing in ([5]) and I is a 
renormalization scale. The quantity Vol represents the volume of the geometrical enclosure that confines the system. 
We have discussed details of the zeta-function formalism for interacting fermion field theories in curved spacetimes in 
a previous publication, Ref. [26j j , and the reader is addressed to that paper for details and additional bibliography. 

The requirement of self-adjointness constrains the form of the boundary conditions and, in turn, the form of the 
function $. The simplest (allowed) choice for fermions is that of bag boundary conditions, 

(l + ih a l )ip\z = , (7) 

where n is a unitary vector normal to the boundary . A general discussion on the boundary conditions can be found 
in Ref. [27j . For the case at hand, starting from (J5J), it is straightforward to show that 



r 



■ fci I + uj^ + a 2 — [i 1 — 2ifiUJ Tl 
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where D = 3 in the case of three spatial dimensiontQ. Integrating over the direction parallel to the boundaries, one 
obtains 

where the first sum is over the solutions of $(fcj_) = and the second sum is over the Matsubara frequencies. In the 
above expression we have Mellin transformed the summand after integration over k\ \ . The last term in the integrand 
can be expressed in terms of the elliptic function , &2{u,v), using the following identity 

n * 

where 

£pA*) = ^1 + 2 £(-1)" cosh(/3/in)e-^^ (10) 
is, in fact, a series representation for i9 2 ,29]. The sum over fc^ in ([5]) is, instead, recognized as the heat-kernel, 

e(i) = ]Te- tfc i . (ii) 

k ± 

The function &(t) encodes the dependence on the boundary conditions and details of the geometry. Putting things 
together we obtain the following expression 

CM = J" dtf- D ^e-^ 2 e(aH)^ 2 t) , (12) 

where, for convenience, we have rescaled t — > a 2 t. The free energy can be obtained by suitable analytical continuation 
of the above expression and its derivative, according to (j6]), once the set-up is precisely specified. In the next sub- 
section, we will consider the case of two parallel layers. 



B. Two-layers Set-up 

In this section we will focus on the case of parallel flat layers that we assume to be located at z = and z = a, 
leaving the field unconstrained in the remaining directions. For the time being, we will specify the boundary conditions 
to be of the form ([7]). We then may proceed by direct computation. The eigenvalue equation may be expressed, by 
decomposing the solutions of the Dirac equation in positive and negative frequency modes and using ([7| . as an implicit 
constraint for the momenta in the z direction, 

(jsin(fc z a) + k z cos(fc z a) = . (13) 

The roots of the above equation are real and positive, consistently with the self-adjointness of the problem. For 
vanishing er, they can be found explicitly, k z = k = ir(n + 1/2) /a with n € N. (A numerical test on the analytic 
approximations for the eigenvalues developed in this and in the next section is reported in Appendix [A| . 

As it should be clear from (I12p . the first necessary step consists in constructing an explicit expression for the function 
O(i). In general, this cannot be done in fully resummed form. However, in several cases approximate expressions can 
be obtained. Partial knowledge of the phase structure can be obtained in the region of high temperature, where it is 
possible to use the small-i asymptotics to get a suitable expansion for the thermodynamic potential. Analogously, the 
case of large a, can also be treated using a similar approximation, since the exponential term in the integrand in f| 12[) 
suppresses the large t portion of the integration range. However, when finite size effects are taken into account, the 
presence of an additional length scale, the inter-layer separation, complicates things, since the relative size of a with 
respect to a must be taken into account. In the following we will construct the analogous of the Ginzburg-Landau-type 
expansion in all relevant regimes and the remaining part of this subsection will be devoted to develop an appropriate 



The case studied in the next section refers to the case of two parallel layers that effectively compactify the direction perpendicular to the 
layers and leave the parallel directions unconstrained. This geometry is clearly non-compact. In this case, we proceed in the usual way 
by closing the cavity with two hyperplanes, separated by a distance L and perpendicular to the boundaries, where we impose periodic 
boundary conditions, and let L to infinity at the end (see, for example, [2£lp. 
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scheme for this. We will present explicit results up to fourth order, but the procedure can be extended to higher 
orders systematically. 

The two regimes of interest are those of large a, relevant to take the infinite volume limit, and that of small a. We 
will begin considering the latter case. The first step of our scheme consists in solving the quantization condition for 
the momenta k± given by Eq. (fT3|) by perturbatively expanding the eigenvalue equation for small a and fixed a. This 
leads to 

/ fj \ 1 fd\ 2 1 „ /a\ 3 1/1 1 ,\ /ct\ 4 5 / 1 4a 2 \ 

where fco &re the eigenvalues corresponding to a — 0. The above expression, ()14[) . can be used in to arrive at an 
appropriate expansion for 0(t). In the regime of small condensate and small a, we obtain 

Q(a 2 t) = @ Q (a 2 t) - 2(aa)G 1 (a 2 t) + 2 {aaf G 2 (a 2 t) - - {adf 9 3 (a 2 i) + | (aaf <d 4 {a 2 t) + ■■■ , (15) 

where we have defined 

9 (a 2 i) = J%(t), 
<di{a 2 t) = tJf (t), 

e 2 {a 2 t) = t 2 ,jeo(t) + ^m(t), 



27 15 

e 4 (a 2 t) = t 4 jeo{t) + zt*jent) - 2t 2 je 1 (t) + — t 2 jr 2 (i) - ztx 2 {t) + -z-tJtr 3 (t), 



and 



jtr a (t) = Y J (ak )- 2a e- ta k o . 

n=l 

Relation (|15[) can then be used to obtain the following expansion for C(s): 



where 



= 


2^(3S 
(4tt) d / 2 


(^Ao - 2 {ad) 


A 


= 4°\s) 




Ai 




i , 


A 2 


= ^\s) 




A 3 


= ^\s) 




A 4 










n 2s-D 

W "r(.) 



3— " 4 



(16) 



2 2 w 2 1 

A 4 = stf' \s) - 2^ 1] (s) + 

with 

dttS _o/2-i+/ e -ta 9 & p>lt (aH)Jtrj(t) . (17) 

For convenience we split the integral (|17j) as the sum of two terms, separating out the thermodynamic contribution, 
Wj J \s) from a temperature independent part, Yj J \s), 

^ J \s) = r} J) (s) + W} J \s) , (18) 

where 

*f )( S ) = / dtt*-^-^e-^ e -^ , (19) 



dttS -I3/2-l+7 e - t a^ 2 -^^^^ (20) 

Notice that both terms depend on the inter-layer separation that in (j2"0|) is entangled with the temperature. In the 
above expressions we have used the following short-hand notation: v = j + 1/2 and u = YlnLi HJLq- Relevant 
representations for the above integrals (|T9"]) and (|2U|) and their derivatives along with the analytical continuation at 
s = are constructed in Appendix [B] for any dimensionality. Composing the results contained in appendix [B] to 
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obtain the thermodynamic potential presents no special complication, aside of being a tedious and long algebraic 
computation. 

A readable expression for the thermodynamic potential may be presented when the various contributions are 
appropriately combined, 

(21) 



^9 



The term comes from the f} J \s) contributions with J = (see appendix IB II part (a)): 



* - Ufi^ t? °" D> {aBT + log (a) s " iD> (<,s) " 



(22) 



The coefficients Un an d t>ti can be computed straightforwardly and for D = 2, 3, 4 are reported in Table (jlj and 

Table (|IT|) in App.0 Notice that the coefficients always occur in finite number and for even spatial dimensionality 
vanish. 

The J 7^ terms can be combined together leading to the following expression ((see appendix IB II part (b))) 



where the first few coefficients are tabulated in Table (IIIII) and Table (IIV[) in Appendix [Cl for D = 2, 3, 4. 
Finally, the temperature-density contribution, can be expressed as (see Appendix IB 21) 



(23) 



%,T = 



2 — a — (- _ 2 / 1 

— — -d <^ woo - 2 (ac) wio + 2 (acr) w 20 + -Q l 
'+i(47t)t >- V 2 



a D + 1 (4 7 r) 

4, , 3 /„ 1. 3. 3, 
--(acr) ( W30 - -u)u + -lo 2 \ + 2^12 

2 . _, 4 / 27 15 

+ 3 ' a<7 '' I W4 ° ~ + ~ + + ywi 3 



)♦■■■>■ 



(24) 



in terms of the functions 



fi„-«£(-i>-(i£ 



n=l 
oo 



I-D/2 



,D-2Z-2J 



cosh(/3/m) x 



i=0 



where the quantity a, is defined as 



©Vl/2)' 



(25) 



(26) 



The results (|22p . (|23p and (|24p can be added up, according to (|2ip . to obtain the final expression for the thermody- 
namic potential. The above method can be systematically extended to higher orders, although computations become 
increasingly more cumbersome. 



C. Expansion for small a /a 

It is easily noticed that in the regime of separation relatively large compared to the value of the condensate, the 
previous expansion is not useful. When the inter-layer separation is larger than the typical value of cr, the expansion 
developed in the previous section does not converge and to study the thermodynamic behavior in the range where 
a I a is small or aa is large is important to develop an alternative expansion. Clearly, this step is necessary also to 
obtain the infinite volume limit. The aim of this section is to modify the procedure to include the case of large a. 
This case can be done essentially in the same way by taking as expansion parameter a /a. After rescaling the proper 
time t —> a~ 2 t in (IT21) . we expand the function 0(f): 



e(t) = e„(t) 





cr" 


9i(t) + 2 


cr" 


2 4 


£7" 


3 2 




2 


.O- 


.O- 


02W - 3 




e 3 (<) + 3 


.a. 



e 4 (t) 



(27) 
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where 



Go (i) = -So(t), 
0i (t) =t£ (t), 

3 (t) = i 3 ^oW + §* a ^i(t) - \<?t£ x {t) + |t^ 2 (t), 

27 



15 



and 



6 4 (t) = i 4 ^o(i) + 3t d ^i(t) - 2a^^i(t) + —rS 2 (t) - 2>a z t£ 2 {t) + —t£ 3 {t) 



Using the above expressions we obtain the following expansion for £(s) 



where 



COO = 


(47T) 15 / 2 ( 


A = 




Ai = 




A 2 = 


^ 0) ( S ) + 


A 3 - 




A 4 = 





An 



Ai +2 



A 2 -^ 



A 



3 T r 



A 4 



(28) 



l fl ^W (s) + ^(i) (s)+ 3^( 2 ) (s) 
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with 



2a 2 SS ( 2 1 \s) +3&g\s) - 3a 2 ^ 2) (s) 
i J) ( S ) = j^y J o to'- D '*- 1+I e-*^(t)3At) 



1^(3), 

2 



(29) 



Here we will present directly the expression for the thermodynamic potential and skip the details of the computation 
of the above integrals that can be done using a procedure similar to the one in Appendix C for the case of small a. 
We express the potential in the following way, 



^9 



(30) 



where +5^q correspond to the sum of the terms with J = 0, *Wi to those with J ^ and to the temperature 
dependent part. Explicitly, + <5%) can be expressed in terms of the function 



rj(s) 



r(s + I-D/2) 
T(s) 



2 x D/2-I-s 
2 



(31) 



where the summation can be written in term of the function 3"{s, p) defined in Eq. (|B7|) . 



ri(s) 



Y{s + I-D/2) fir 

r£) 

-*(a + I-D/2, ^ 

V 7T 



2X D/2-I-s 



-,2s+2I-D 



^ S + I-D/2, 



2aa\ 



(32) 



Proceeding as described in Appendix IB 11 we get 



r 7 (s) 



r(a) Va 2 



2X D/2-I-s 



*-T(s + I-D/2 -1/2) — 



l-2s-2/+_D 



+ 2tt 



i+I-D/2 



X 

V 7T / 



-3=1 

Ka+i-o/2-1/2 (2qaa)] } 



-D/2-1/2 



2 s+7 - D/2+1/2 ^ +/ ^/ 2 ^ 1/2 (4gaa) 



(33) 
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The terms and <5%) collect the contributions coming from the first and second terms in curly brackets in the 

expression above, respectively. It takes straightforward steps to arrive at the following formulae: 



(i) 



2 — 



1 



(r {aa 



D+l 



(47T) 



2r7{aa) D + 2rJ (aa) 



D-l 



o r 3 ( acr ) 



D-2 



r+(acr) 



D + l 



2r+(aa) D + 2r+(aa) 



D-l 



3 J ' ' 3 
where we have adopted the short-hand notation 

D + l 



= — + r+ + O(s) , for s -> 



(34) 



(35) 



Notice that the coefficients r7 vanish for even number of spatial dimensions (even D), while only the first (D + l)/2 
are non- vanishing for odd number of spatial dimensions (odd D). We have defined a e = exp (j e /2) with j e being the 
Euler constant. The contribution comes from second term of (|B11[) and it can be expressed as 



2^ 1 



_ -i 1 f 4 Z 

(47T) 2 M *. O O 



where 



-^H^^^-^ [2 I - D l 2+1 / 2 K I _o+ 1 (Aqaa) - Kj_o^ (2qaa) 

V 0=1 



(36) 



(37) 



The J terms can be combined together leading to the following expression 



2 — 



1 



(47T) 13 / 2 a D+1 

-(f) Va 



m + 2 ( - 



1 2 

-a Tn — T 2 i — Ti2 



-T 2 i ~ T31 + a Tl2 — -r 22 — -Ti3 



where tjj is constructed from the analytical continuation to s = of 



u 



i D+2J 

^2J 



r(J-£>/2 + s) 

r( s ) 



-2 i " 2 



D/2-I- 



and its derivative, giving 



s->0 



r /J = lim ^jlogr 2 + — sr u 



(38) 



(39) 



(40) 



A regular integral representation for the above series can be obtained by use of the Abel-Plana summation formula 
(see appendix B for details): 



tij 



^(f) ' & D/2 - I (x + -X-log{Z 2 &))+a D ^-) 2 \ IJ 

(z + /3o + ^/3i-^/3 log(ftr) 2 ) 



2(1+ J) 



2- 2 n(1+-D-2,7-27)/2 

-a a 



where & = £?>(0) with 



&>(£) 



k z 1 



and 



^13 = I 



\ + w) " ^^V) (x- log (r a ^( W )) + x+) - (z -*)] ^^J^ 
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where the constants x±: z ± an d A) an d Pi are defined in (ID9|) . (ID13I) . (|D14[) and (|D15|) . 

The above integrals can be evaluated analytically by dividing the integration range into two intervals, [0, 1) U [1, oo] 
and by expanding the exponential appropriately. Alternatively, use of numerical approximation does not require any 
effort. 

Finally, the temperature-density contribution, can be expressed as 



where 



%,T = 



n D — l 

2~ 



±Z1 J 

> D + 1 (47T)T I 



4 


cr" 


3 ( 


~3 


-O- 




2 


a~ 




+ 3 


.a. 





1 



13730 



"<7~ 






<7~ 


2 f 






wio - 






1 ^20 "I 




-Q _ 






.a. 







2 i J 

-a ron + 2^21 + 2^12 



CC740 — 2a cc72i + 3tx73i — 3a ra7i2 



27 



-G7 2 2 



15 



-tUi3 



(41) 



2I+2J- 

a ujij . 



(42) 



Combining (IMl) . (|3^|) . ([35]) and (|4ip according to ([2^1) and (JB]) gives the thermodynamic potential. 



D. Infinite volume limit 



In the present section we will explicitly take the infinite volume limit. As we shall see, aside for the cancellation of 
spurious divergences, the computations that follow will provide a non-trivial check of the general results. 

To explicitly take the infinite volume limit, rather than the above integral representation we can obtain a simpler 
expression by using the binomial theorem. Choosing aa <C A € N, we have 



a D ^ T{I- D/2 + s) 

JlJ = — rr^ X 



(D/2-I-s + k- 1 




- 1\ / 7T \ 2fe 

1 b) CH(2J-2 fc ;l/2) 



2fe-2J 



(43) 



Starting from the above expression, simple power counting arguments are sufficient to show that the infinite volume 
limit can be taken with no problem. The last two terms involve a finite summation over j [y = j + 1 /2) and an infinite 
summation over k. Both of them do not produce any diverging behavior, as trivial algebra shows. The first term 
is also regular, since the zeta function occurs with an even argument. Both summations over k are also converging 
in the infinite-volume 'regime', since aa is large. However, using the above representation directly in (|38[) and (|39p 
requires some care when the limit a — > oo is taken since the leading behavior of the above expression for large a goes 
as 3?h ~ a 2J and not all terms in (|38|) vanish in this limit. Precisely, in the second line of (|38|) the terms proportional 
to a 2 Ti2 and — 5T13/2 individually produce linearly diverging contributions as a approaches infinity. These terms, 
when combined, cancel out leaving a regular expression: 



r(s) 



C(4;l/2) C(6;l/2) 



n=0 



n=0 



n=0 



n=0 



0(l/a) 



= + O(l/a) 



(44) 



Although the representation (|43|) is sufficient to prove the regularity of the infinite volume limit, for large, but finite, 
a the integral expression given above is handier for numerical evaluations. 
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Letting a — > oo, we have already shown (see (1441) ) that the second line in ^i, Eq. (l38l) . vanishes. The first and 
fourth terms in the first line of also trivially go to zero, while the second and fourth asymptote to constants 
that when combined vanish. The term S^o also vanish in the same limit. The contribution ^.t, is more delicate 
to analyze. All terms in the first line in (1411) asymptote to zero in the infinite volume limit. In the second line of 
the same expression the first and second terms vanish in the a — > oo limit, while the second and third asymptote to 
a-indcpendent quantities. In the third line of (|41j) . the first, second, third and fifth terms go to zero. The forth and 
sixth terms, while individually diverging as ~ a, when combined, vanish in the a — > oo limit (the reader may check 
this by expanding the quantity — 3zui2 + Vhwnj2 for large a and perform the sums over i exactly). All terms that 
asymptote to a-independent quantities in the infinite volume limit, vanish when the zero temperature limit is also 
taken. Therefore, for T = 0, only the first term coming from the contribution survives in the a — > oo limit, 

giving 

Hm % W = 2 ^ D+1 * D+1 T ( 1 - °p^\ . (45) 

The above result is seen to reproduce the correct infinite volume zero temperature limit (see Ref. [20T ] formula (31) 
or Ref. (30jV The terms in ^.t, that do not vanish in the infinite volume limit, reproduce the correct temperature 
dependence as it can be checked, with some work, by comparing with the flat space limit of the results of Ref. (26j 
(see also 0]). 



III. REMARKS 



In the present section, we wish to make several remarks useful for the subsequent numerical evaluation, as well as 
perform some other checks. 

Vanishing Condensates. The limit of vanishing condensates, a = 0, at zero temperature and chemical potential 
provides one non-trivial check of the result. We will consider the case of D = 3 in order to compare with the correct 
zero temperature free energy for masslcss free fcrmions confined between two parallel layers at distance a given in 
Ref. [HI 

7tt 2 

W c = ■ 7 . (46) 

2880a 4 V ' 

This result should be reproduced by our computation in both regimes of small and large separation. In the case of 
small separation, the reader may easily check that the cr-independent contribution comes from the n = term of the 
first summation in (|22|) . This (see Table ((XJ) for the numerical coefficient) reproduces (|46|) . For large separation, the 
leading a-independent term may be computed by expanding (|36[) and by performing exact resummation over q. This 
also reproduces the above result (|46|) . Finally, let us notice that the same result can be obtained directly setting a = 
in (|12[) . In the limit of vanishing temperature and chemical potential, a simple computation gives 

-^0^0.0.-^(^-1). 

where the function 5f is introduced in (|B7[) . The above expression can be inserted in ©, and reproduces the correct 
result. 

Summations. The results obtained in the preceding sections contain series over certain combinations of Bessel 
functions. Here we will examine the convergence of these series and show how to implement their numerical evalua- 
tion. We will illustrate the procedure for the case D = 3, extensions to any dimensionality being identical. Expression 
([57)) can be written in terms of the quantity 



oo oo 

= J2 & q = Yl I 1 ' 2 ( 22_/ ^-2 (4<?a<7) ~ Kx-2 (2qaa)) 



(47) 
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The above sum may be evaluated in a way that is numerically efficient by separating out the small-g part, 



oo oo 

= £^JL + E {^i-^l) > (48) 



q=l 9=9+1 



with integer q ~ (eta) 1 and the underline signifying expansion for large argument. The first term in the above 
expression can be resummed exactly to any desired order. Each term may contain spurious divergences that can be 
handled by multiplying the summand by q u before resumming over q and then by analytically continuing to u = 0. 
The second term constitute the reminder of the summation and it is expected to be small due to the exponential 
decay for large arguments of the Bessel functions. Computation of the reminder can be handled in a numerically 
efficient way by using the Abel-Plana summation formula (see, for instance, |32j), 



v i />oo pC 

V Jf{q) = -Jt?{a) + / Jf(z)dz + i \ 
„_„ 2 J a J Q 



J^(a + ly) - J^la - ly) , 

^73! dy 



The reminder is found to be negligible in all cases analyzed. 

Numerical Evaluation of the r-functions. Apart from the integral representation for the functions tjj(s) given 
in the preceeding section, it is possible to evaluate these efficiently in a fully numerical way. A convenient approach is 
to fix the values of D, I and J first, then expand r/j(s) and its derivative around s = 0, and then sum the resulting 
expressions numerically. Numerical summations can be performed in a number of ways, here we adopt a sampling 
technique based on computing the summations up to a cut-off, sample the remaining terms, and approximate using 
a polynomial interpolation. The reader may easily check that for I + J > 2 this numerical approach works very 
efficiently. However, something different is necessary to deal with the case I = J = 1. In this case, we adapt a 
technique similar to that used in Ref. [33| to extend the Chowla-Selberg formula for the Epstein zeta function to 
higher dimensionality. The trick is to write the second term in the summation occurring in (1391) in terms of its Mellin 
transform, discretize the o direction and integrate over t at each point of the a-grid (the exponential dependence on 
the condensate allows one to interchange the summation and integration operations - see Ref. [33[ for details about 
this point in general). The resulting expression can be easily used to construct an interpolated form for tjj. 

Convergence. An analogous scheme may be used to compute the function uji t j occurring in the temperature-density 
contribution, however, due to the presence of a chemical potential, the analysis of the convergence of the n summation 



is more delicate. Simple steps show that the series converges exponentially for fi < y + & 2 setting a limit of 

validity for the series representation ([25]) . A rather nice consequence of the above inequality is that by decreasing the 
inter-layer separation the convergence improves and allows for larger values of [i than in the corresponding infinite 
volume limit. This is rather convenient for any numerical evaluation, the reason being that for values of a and fi 
satisfying the above inequality the sum over n can be approximated by keeping only the first few terms, and any 
numerical computation does not require any effort. For small [i it is possible to expand the summand appropriately 
and resum each term of the expansion by means of analytical continuation techniques. 

Inter-layer Separation. It is interesting to sec that the inter-layer separation, a, appears in the expression for 
the thermodynamic potential in a form identical to that of the Matsubara frequencies (|26p with an effective tem- 
perature T a = l/(2a). The presence of a chemical potential and higher order terms in the condensate makes the 
combination between geometrical and thermodynamical effects non-trivial, making the scrutiny of the finite-size - 
finite temperature analogy not as transparent as in the case of geometries with topologically non-trivial circular 
directions [22]. 

Scaling. Letting [mass] — > K[mass], where [mass] is any mass parameter (a, /i, g -1 / 2 , /3 _1 , a -1 and it 
takes simple steps to show that the potential rescales as °i/ — » K D+1 ^f , allowing one to fix one of the parameters. 

Large condensates and high temperatures. The high temperature approximation for the potential in the regime 
of large temperature can be obtained with less work than we did above, by using the method of Refs. [H], |35| which 
consists in rescaling the proper time in the expression for the zeta function (fT2"]) , t — > t/3 2 /(27r) 2 , and then uses 
the small-t (Schwinger-De Witt) expansion for the heat-kernel. A similar trick can be used to analyze the large 
condensate behavior, as it is easy to understand from the general formula of the zeta function. Rather than reporting 
this computation for the two-layer set-up, which is essentially a repetition of the previous calculation, we will consider 
both cases in some detail in Appendix [El for more general boundary geometries. 
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IV. RESULTS 



The discussion presented in the previous sub-sections, along with the details reported in the appendices, should 
provide sufficient information on how to proceed with the numerical computation. A technical complication comes 
from having to use different representations in different regimes, but, apart from this, the procedure is straightfor- 
ward. We will present the result of such an evaluation for the 3-dimensional case and 2-dimensional boundaries. 
Generalization to different dimensionality can be performed similarly. 

Zero chemical potential. In the first part of this section we will consider the case of vanishing chemical potential. 
We have already explained the basics of our method in the previous section. One thing that deserves further discussion 
is the way we perform the summations over the Matsubara frequencies and over the eigenvalues o^. Clearly, both 
sums can be performed numerically using standard numerical techniques. However, a more efficient way to proceed is 
to expand the Bessel functions in the region of interest and then exactly resum over the Matsubara frequencies each 
term of the resulting series. This provides an analytical expression that will need to be summed over the eigenvalues 
cti . This second sum is then performed numerically. The full numerical method and partially resummed one are found 
to give consistent results. 

In Fig. [I] we illustrate the result for the thermodynamic potential ^Ov (normalized by subtracting the value of the 
free energy for o — 0) for sample values of the parameters. Notice that when the inter-layer separation is decreased, the 
effective temperature of the system increases and the minima of the potential shifts towards smaller values. Another 
feature that can be noted is that the potential 'flattens' when the value of a increases. The third prominent feature is 
that for combinations of distance and temperature it is possible to arrange the system such that the effective potential 
has two minima (one vanishing and one non- vanishing) with equal free energy, suggesting that in certain regions of 
the phase diagram, possibly implying the existence of regions of broken/restored phase separated by domain walls. 
In this case, the condensate would not be constant but spatially varying, thus leading to an inhomogeneous phase. 
Whether, even at vanishing density, such a phase could become energetically favored owing to finite size effect cannot 
be proved at this stage, but we will return on this point in the future. 




FIG. 1. The figure displays the thermodynamic potential (normalized according to ^| s =o = 0) for vanishing chemical 
potential as a function of the condensate o for fixed inter-layer separation and varying temperature. The values of the parameters 
are fixed to g — 10 2 , l~ x = 30. The value of the separation a is indicated in the panels, while the temperature takes the following 
values: = 7. (top curve - left panel), = 1.5 (bottom curve - left panel); = 4.5 (top curve - central panel), = .7 
(bottom curve - central panel); /3 _1 = 2.2 (top curve - right panel), = 0.5 (bottom curve - right panel) 

In Fig. [2] we illustrate the a-dependence of the potential for fixed values of the temperature. The dashed red line 
represents the critical curve at which a phase transition occurs (for fixed temperature) owing to finite size effects. 

In all figures, increasing the inter-layer separation or decreasing the temperature (for fixed small separation) tends 
to flatten the potential and eventually change the order of the transition from first to second. 

Finite density. Switching on the chemical potential produces some complications related to the fact that its 
dependence modifies the properties of convergence of the expression (|25p . However, as we have stressed previously, 
the incorporation of effects of finite size compensates that of a non- vanishing chemical potential allowing for values of 
fi larger than the corresponding ones in the infinite volume case. The various regimes that can be considered may be 
treated more efficiently using different numerical schemes. For instance, assuming to be in a range of the parameters 
for which the convergence of [[25]) is guaranteed, if n/T <C 1, it is convenient to expand the summands in (f25j) for small 
chemical potentials, proceed by exactly resumming over i and then sum over the Matsubara frequencies by means of 
numerical approximation. When the chemical potential is large compared to the temperature, it is more convenient 
to perform both summations numerically. Here, we always follow this second approach and use the first only for 
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FIG. 2. The figure illustrates the dependence on the inter-layer separation of the thermodynamic potential for fixed temperature. 
The values of /3 are set as indicated in the panels, while the separation a takes the following values in all panels: a = 5 (top 
curve), a = 20 (bottom curve). The red dashed curves represent the critical curves: (left panel) a — 10.4, (central panel) 
a — 13.2; (right panel) a = 15.12. The other parameters are set as in the previous figure. 

check. The double numerical sum is performed over a grid in the a direction. We perform exactly the sum over the 
two indices up to a cut-off value, sample the additional terms and approximate these using a polynomial fitting. The 
values of the cut-off are varied until the desired accuracy is reached. Results are shown in Fig. ([3]) for sample values 
of the chemical potential and of the inter-layer separation. 




0.0 0.5 1.0 - 1.5 2.0 2.5 0.0 0.5 1.0 1.5 -2.0 2.5 3.0 3.5 4.0 1 2 - 3 4 5 

a a a 

FIG. 3. The figure illustrates the thermodynamical potential for vanishing vs non-vanishing fixed a and varying /?. The 
inverse temperature is varied between /3 _1 = 0.6 and /3 _1 = 3 for the left panel, /3 _1 = 1 and = 4.5 for the central 
panel, = 3. and /3 _1 = 7.5 for the right panel. Values of fj, and a are indicated in the figures, coupling constant and 
renormalization scale are fixed as in the previous figures. 

The critical temperature, T cr i tl can be computed by solving the equation 

= lim d 2 W/8a 2 = (49) 

er— >0 

by numerical approximation. We first expand to second order in a the thermodynamic potential, and, after computing 
the second derivative of the resulting expression, we take the limit a — > 0. The infinite summations over the Matsubara 
frequencies and over the eigenvalues ctj involved are computed using the same numerical sampling technique we have 
used before. The sums are evaluated after discretizing the /3-direction and for fixed values of the chemical potential 
and of the inter-layer separation. The result is then fitted using a polynomial interpolating function and the zero of 
fa,n{P) = is calculated numerically by using the interpolating function. Alternatively, the critical temperature can 
be computed directly by inspecting how the minima of the potential changes when the temperature and chemical 
potential are varied and by finding for which value of the temperature for fixed /j, and a the non-zero minima of 
the potential vanishes. This method is accurate when the transition is first order while it is less accurate when the 
transition is second order. Results are presented in the right panel of Fig. 

We close the section by computing the value that the condensate attains at the minima of the potential, a m i n , vs 
the temperature (for fixed separation) and vs the separation (for fixed temperature) both at zero chemical potential, 
and the minima are computed using a minimization routine based on the Newton's method. The left panel of Fig. ((4]) 
below shows the value of <J m in obtained for a sample choice of parameters. 

A precise construction of the three dimensional phase diagram T-/i-a is left for future work as it requires a more 
substantial numerical effort (particularly for the precise determination of the critical points and of the order of the 
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FIG. 4. The left panel illustrates the behavior of a m in vs the temperature for several values of the inter-layer separation (values 
are indicated in the figure). In the right panel we display the critical temperature, T cr n vs [i for several values of the inter-layer 
separation. The filled dots are directly calculated as explained in the text, while the empty dots are computed by fitting the 
results in the region around the 'knee' and by extrapolating the resulting curve. The renormalization scale and the coupling 
constant are set as in the other figures. 

transitions). However, the results displayed in Figs. (j2]-{4j) are sufficient to illustrate the main features that the presence 
of the layers produces. As expected, reducing the inter-layer separation induces a shift in the thermodynamic potential 
indicating an increase in the effective temperature. This is also observed in the tendency of the critical temperature 
to decrease when the separation decreases, indicating that for fixed (3 a decrease of the inter-layer separation tends 
to bring the system towards a phase of restored chiral symmetry. An additional interesting feature is noticed looking 
at the (3 dependence of the potential for fixed a and varying [i (see Fig. ([3])), showing that a simultaneous decrease in 
temperature and separation tends to cause a flattening of the potential and eventually a change in the order of the 
transition from a second order one for large a towards a first order one for small values of a, and producing possible 
shifts in the critical points. 



V. CONCLUSIONS 

It is not possible, in general, to predict how the finite size of a system may affect its thermodynamical behavior. This 
depends on the geometrical and topological properties of the system, on the presence of external fields, on the specific 
type of interactions, etc. In some cases, symmetries may help to relate high and low temperature-regimes, or (modular 
symmetries) certain topologically non-trivial configurations. Spectral geometric techniques can also be used to make 
more generic statements in specific regimes {e.g. at high temperature). However, a case-by-case analysis seems to 
be necessary to make detailed predictions and understand precisely the combined geometrical and thermodynamical 
effects. 

Here, we have analyzed the impact of finite size effects on the thermodynamics of a system of interacting fermions 
that we have modeled using an effective field theory with a four-fermion interaction term in D dimensions, and 
looked in detail at the case of a confining two-layer system in three dimensions. We have developed a method based 
on zeta-function regularization techniques and used the large-iV and mean-field approximations to construct the 
effective thermodynamic potential. We have developed all necessary expansions in different regimes of temperature 
and separation, constructed all relevant integral and series representations for the various expressions involved and 
thoroughly discussed the analytic structure of the potential. While being particularly interested in the case when the 
'tree-level' four-point interaction term and quantum corrections produced by the deformation in the vacuum due to the 
presence of the boundaries both generate a sizeable force between the boundaries, we have also obtained an expression 
valid for large separations that is important when taking the infinite volume limit. We have performed several checks 
of the individual contributions showing non-trivial cancellation of divergences and reproducing various known results 
in several limits (infinite volume limit, zero condensate, finite temperature). The method we have discussed in this 
paper can be efficiently implemented numerically and, after clarifying several points, we have performed a numerical 
computation of the thermodynamic potential, of the critical temperature, and the condensate in the broken phase for 
the D = 3 dimensional two-layer set-up. 

The most important effect of finite size is the tendency to shift the critical points and eventually change the 
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order of the transitions. We have seen that the critical temperature and the 'flatness' of the potential depend in a 
non-trivial way on the size of the system (the inter-layer separation in our case) and on the other thermodynamical 
quantities. However, in order to completely describe the modifications that finite size effects may produce to the phase 
diagram of a strongly interacting fermionic system, it is important to relax at least two conditions. The first one is the 
approximation of spatially constant condensates, while the second one is related to the boundary conditions. Relaxing 
the first condition generates several new terms in the partition function and solution of the gap equation in this case 
requires proper functional minimization and related numerical treatment. Regarding the boundary conditions, here 
we have considered the simplest possible choice allowed. However, more general boundary conditions may be used 
[27j . How the phase diagram depends on the boundary conditions remains an open (and in our opinion interesting) 
question and it is clearly very important also in the inhomogeneous case. When the boundaries are disconnected, as 
for instance in the two-layer case analyzed here, different boundary conditions may, in fact, be imposed at each layer 
independently (introducing in the effective potential a dependence on some additional parameter, e.g. a chiral angle). 
While these changes are not expected to modify the situation at high temperature (see for instance appendix [EJ some 
modifications may be expected at low temperature. Especially when inhomogeneous phases become energetically 
favored, how the ground state is affected by the boundary conditions is certainly non-trivial, and whether by means 
of an appropriate choice for the boundary conditions it is possible to engineer the shape of the ground state is also 
an amusing related question. 



massless phase 
F c ~ a" 4 



massive phase 
F c ~ a- 4 e" a5 " 







FIG. 5. VEV of the fermion condensate and phases of the Casimir force. The figure shows how the VEV of the condensate 
depends on the inter-layer separation in the region around the transition for fixed temperature (set to = 1.5.). The other 
parameters have been set as in the central panel of Fig. ([2]). 



A physical situation where the present discussion may have some relevance is in the context of the Casimir effect. 
Here, we have basically considered a Casimir-type set-up with the action augmented by a quartic interaction term. 
An interesting system of this sort is the bilayer graphene with interlayer pairing whose properties have been rece ntly 
discussed in Ref. (36|. (The literature on the fermion Casimir effect is vast. The reader may consult Refs. [37H40j 
for some recent articles where the Casimir effect for non-interacting fermions has been discussed.) While in the free 
field case, no phase transition/symmetry breaking occurs, adding an interaction term may change things in a more 
dramatic way. Let us choose the coupling constant g and the temperature in such a way that for large separation a 
chiral symmetry is broken (as in Fig. Starting from this configuration, let us gradually decrease the inter-layer 

separation producing a gradual shift in the minima of the condensate. Once the distance reaches a critical value at 
which chiral symmetry is restored a phase transition occurs: the true minimum of the potential a m i n jumps to a 
zero value (see Fig. [SI where the transition is a first order one). The transition is from a massive phase (with a mass 
proportional to the condensate) to a massless phase in the restored chiral symmetry region. In the massive phase, the 
fermion contribution to the Casimir force is expected to be suppressed (with the degree of suppression depending on 
the precise value of aa m i n ). On the other hand, in the massless phase, the fermion contribution to the Casimir force 
takes the familiar a~ 4 behavior. The transition between these two phases of the Casimir force (mass-suppressed and 
massless) is related to the spontaneous breaking/restoration of chiral symmetry and it is a phase transition. Fig. ([5]) 
summarizes the above arguments. 

It seems interesting to us to pursue these ideas further and we will return to these issues in the near future [23j . 
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Appendix A: Numerical test on the analytical approximation for the eigenvalues 



In order to test the approximations we have developed for the eigenvalues in Sec. Ill Bl and Sec. Ill CI we have 
numerically computed the eigenvalues and compared the resulting values with the analytic approximations over a 
wide range of parameters for k± up to 10 5 . Fig. [5] and Fig. [7] show some sample plots of this test for small and large 
values of a and a. The intersection between the vertical red lines and the horizontal axis indicate the eigenvalues 
computed by numerical approximation using a root-finding routine based on the Newton method. The blue dots are 
the eigenvalues obtained by analytical approximation up to fourth order in the expansion parameter. 
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FIG. 6. Numerical check on the eigenvalue approximation. Plots refer to a = 0.1 and a = 0.1. 
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FIG. 7. Numerical check on the eigenvalue approximation. Plots refer to a = 1 and a — 10. 



Appendix B: Integrals 



In this section we will provide some explicit formulas for the integrals used in the body of the paper (some for- 
mulas and definitions given before will be repeated here for the convenience of the reader). The starting point is the 
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following expression 



where 



Jls — D poo 

*f)(„) = ^_/ dtt s - D / 2 - 1+I e- ta2 * 2 J?^(a 2 t)Jtrj(t) , (Bl) 
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= ( B2 ) 



and 
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n=l 

We have used the short-hand definition v = j + 1/2 and = X^o- We wm S P^ ^ ne integral (IB1|) as 

^ J ) ( s ) = y« ( s ) + J) ( s ) , (B4) 

where 
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1. Computation of "V^' r \s 

The first integral can be computed straightforwardly leading to 
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Let us now consider the two cases J = (a) and ,7/0 (b) in turn, 
(a) For J = we have 
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where we have expressed the sums in terms of the following function 



y( S ,p):=J2(n 2 + p 2 y S ■ (B7) 

n=l 

The cases a = and a ^ present some slight computational differences. For a = (p = 0), S^(s,p) is nothing but 
the Riemann zeta function, .3^(s,0) = Cr(2s): 

r/ 0) ( S ) - a 2s - j r(/ ~ ^ /2 + 5) ^- 2s - 2J+g (2 2s + 2/ -^ - 1) Cfl (2s + 27 - D) . (B8) 
1 (sj 
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For any even spatial dimensionality (any integer and even D), the zeta function is regular and the only poles occur in 
the Gamma functions for s = (denominator) and for s = and 1 — 0,1, • • • , -j (numerator) compensating each 
other and leaving a well behaved expression. For any odd spatial dimensionality, any integer and odd D, the zeta 
function has a pole occurring for s = when / = (D + l)/2. This pole is compensated by that of the gamma function 
in the denominator, and the expression is perfectly regular. 

When both a and a are small (that is when p is small), we may expand the function .y(s,p) as a series of Riemann 
zeta functions, 

y(s,p)=^{~iy{" + ^Ap 2 ^ R {2 S + 2q) . (B9) 

Using this representation we get 



r<f» fQ s _ „2 S - D m-D/2 + s) _-2 S -2i+D VV,™ fs + I-D/2 + q-l 
1 1 j T{s) Q 1} { s + I-D/2-1 



x ^2s+2i-D+2 q _ 1) ^ 29 Cfl ( 2s + 2I-D + 2q) 



that recovers ()B8|) in the appropriate limit. As before, for D even and integer, the poles of the Gamma functions 
compensate each other leaving a regular expression for s — 0. For D odd and integer, the poles in the summand, at 
s = and for / = (D + l)/2 — q e N, are canceled by the pole of the Gamma function in the denominator, leaving a 
regular expression. 

A representation valid also for larger values of aa can be constructed by using the Chowla-Selberg representation for 
the series, Ref. |4lj . 

^(s ! P) = 2~ + ~2 f(s) 9 T(7j P 2-^ Q K s-i/2 {2nqp) . (BIO) 

Using the above expression, we find 
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9=1 

-K.+1-D/2-1/2 (2gaff)] } . (Bll) 

One may notice that the first term in square brackets in (jBlOj) cancel out in the above expression. For D odd and 
integer, the poles of the gamma function in in square brackets are compensated by the pole of the gamma function in 
the denominator, leaving a regular expression. For D even integer, the expression vanishes for s = 0. 
(b) For J ^ 0, we have 

*f>( S ) = a 2s- D ni- 0/2 + 8) ^ y ^ 2J (fl2 _ 2 j j y.- I+D/ 2 

i r(a) v J 



For vanishing cr we have 



~l/( J )( a \ _ „2i-D ^ D/2 + s) s—^ 2s-2I+D-2J -2S-2I+D-2J 

= T(I-D/2 + s) __ 2s _ 2I+D _ 2J 
where 

OO 

C ff (w;6) = ^(fc + 6)- tl (B13) 



r . , Cff(2s + 2J-£> + 2J;l/2) , (B12) 

r(s) 



fc=0 



is the Hurwitz zeta function. Since the only (simple) pole of the Hurwitz zeta function (|B13I) occur at u = 1, for 
any D even integer the Hurwitz zeta function in (|B12I) is regular for any positive integers / and J, while the poles 
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coming from the gamma function in the numerator for I = 0, • • ■ , D /2, are compensated by the poles coming from 
the gamma function in the denominator. For D odd integer, the poles occurring in the Hurwitz zeta function in (|B12|) 
are compensated by those from the gamma function in the denominator. As in the previous cases, the final expression 
is regular. 

When aa is non-vanishing and small, use of the binomial theorem leads to 



y{J) ^ = a 2s-D -D/2 + S) ^ 2 I+D-2s-2J 



DO 

q=0 



T(s) 



(B14) 



2q 



The regularity of the above expression for s = and for any non zero integer p can be shown as in the previous cases. 
It is also possible to give an integral representation for "f^" 1 ' (s) using, for instance, the Abel-Plana summation formula. 
However, we won't present explicit formulas for this case here. 



2. Computation of W^ J) (s) 

Wj J \s) is easier to handle and simple steps lead to 

I+s-D/2 



T(s) ^ 

v 1 n—l 

^-2J \ ^ ,,-2J / _ 
7T y V (7T 



cosh(/3/in) x 



2 2 , 2-2\-k(I+s-D/2) 



(B15) 



K I+s _ D/2 (n^(^ + aWf 2 ^) 



The presence of the Bessel functions in the above expression suppresses, for large values of their argument, the 
summations and make any numerical evaluation straightforward. The regularity is immediate to show and the above 
expression (but not its derivative) vanishes in the limit of s — for any D, I and J. 



Appendix C: Tables of the coefficients Uq D \ Uq D \ v\ D \ 
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TABLE II. Coefficients t& 
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TABLE III. Coefficients u, 
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TABLE IV. Coefficients v, 



Appendix D: Integral representation of the .^./-functions. 



Carrying out the analytical continuation of the functions S?u occurring in the regime of large inter-layer separation 
requires an approach slightly different from that we have used in the other cases. The starting point here is the 
function (we suppress the indices / and J for notational convenience and write ST = 



-2.1 



T 2J 



T(s) 



n=0 



D/2-I-s 



(Dl) 



The case of positive J £ N can be dealt with by using some recursive relations starting from the standard Chowla- 
Selberg representation (see 42]). Unfortunately, here we need to consider the case of — J € N and this complicates 
things slightly. Using the Abel-Plana formula, we can express r as 



where 



ST = 



T 



(o) 



ji(O) = 



a 2J T(I - D/2 + s) 
Z2J 



r( s ) 



2 2J-1 ggD/2-I-s 
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T (2) = « 2J r(J - D/2 + s) 
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r( s ) 
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where we have used the short-hand notation 
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It is easy to get 
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T (o) <?°/ 2 - I [x- + s( x+ -x-loz 



+ 0(s 2 ) , 
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where we have defined 



r(J - D/2 + s) 

r( s ) 



X- + S X+ + 0(s 2 ) 



(D9) 



Notice that \- — for D odd. 

In general, the contribution can be expressed in terms of hypergeometric functions. However, a representation 
more useful for our purposes can be constructed in terms of incomplete beta functions, 



p a (x,y) = / t*- l {l-ty-Ht , 



and their regularized form, f3^ eg (x , y) , 



P? g {x,v)=0 a {x,y)/0{x,y) 



Some computations lead to 



= a 2(/+.7+ s)r(/ _ D/2 + s) (1+D _ 2J _ 2J _ 2s)/2 
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D + l D 
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r r D + l D r 
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where the incomplete gamma function is defined as It takes simple steps to show that the quantity 



T(I - D/2 + s) 



f3[D/2-2(I-l + s),I + J + s- 



D + l 



= z-+ z+s + 0(s 2 ), 



is regular for s — >• 0, while the rest of (|D12|) is regular by construction. Expanding for s — > gives 

a 2(I+J) 
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2tt 
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where we have defined 



d 
ds 



Finally we deal with T^: 



(-i) 1+D - 2 *- 21 - 2 ' ^ re 3 s . 
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(D12) 



(D13) 
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where 



J = i 



-2J 



+ iz) & D / 2 - z (tz) - (z -z) 



dz 



- 1 



(~ + «) ^ ^ >D/2 - I ^z) log 9{%z) - (z -z) I 



Analytical approximations may be obtained by splitting the integration range into [0, 1) U [1, oo] and by appropriately 
expanding the exponential in the region of small and large z. Both integrals J? and J! can be computed with no 
effort by means of numerical approximation. 
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Appendix E: High temperature and large condensate approximations for general smooth boundaries 

In the present section we wish give some more details on the high-temperature and large-condensate expansions. 
In these cases, it is possible to give simpler and more direct derivations as well as considering slightly more general 
geometries, for the cases where explicit knowledge for the heat-kernel expansion is available. Although it is possible 
to include some singular spaces (for instance conifold geometries, see [43j), here, for simplicity, we will limit our 
discussion to the case of co-dimension one, non-singular boundaries. The starting point, i.e. Eqts. ([S])-®, remains 
valid, and the zeta function takes the following general form 

« s ) = E E (k - + A + ~ s > ( E1 ) 

n X 

where A is used to indicate the eigenvalues of the Laplacian with boundary conditions appropriately imposed. As we 
did before, one may easily recast the above expression in factorized form 

C(s) = ——r=——— f°° dtf-^e-^e (t) #^ (t) , (E2) 

where 0(f) = £ A e~ tA . 



1. High Temperature 



The first regime we wish to address is the high temperature one, where we expect the system to be in a region of 
unbroken symmetry. In order to study the high-T regime, we may rescale t — > -^-st, in ()E2|) and arrive at 
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At high-temperature (/3 -C 1), we may use the asymptotic expansion of the hcat-kcrncl 
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where the coefficients 9^ are determined by the intrinsic and extrinsic geometry of the spatial section of the background 
manifold as well as by the boundary conditions. The high temperature expansion, skipping lengthy computational 
steps, can be cast as follows 
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The coefficients Vk calculable numbers (for D — 3, the first few are: vq — 3/4, V\/ 2 — 4y / 7r/3, Vi = —1, v 3 / 2 = 

2y / 7t, v 2 = 0, • • • ), and we have defined §k = VoP 1 ^.. 

The above expression is valid for any smooth boundary and for any self-adjoint preserving boundary conditions, 
encoded in the integrated heat-kernel coefficients. Notice that the chemical potential, in the high temperature ex- 
pansion, shows up only at third and higher orders and it multiplies volume (at third order) and boundary (at fourth 
order) terms, but not the condensate (the same is true for bosons 22]). Therefore, up to fourth order in the high 
temperature approximation, it does not affect the position of the minima of the potential. 

The explicit form of the heat-kernel coefficients is not known in general. However, for the present analysis we only 
need information regarding the first few coefficients and this can be retrieved from general formula? (see, for example, 
El])- In fact, the relevant behavior of the potential can be obtained by looking at the sign of its derivative 
with respect to a, e = Sign^||-. For the case of bag boundary conditions, the first few coefficients are known and 
sufficient to show that e > for any positive ct, implying that the potential is a monotonically increasing function 
of the condensate with an absolute minima at a = 0. Therefore at high temperature and without boundaries chiral 
symmetry is always restored at leading order in the high temperature approximation, as one may have expected. 

For more general boundary conditions (for example, for chiral boundary conditions) the first few coefficients have 
been computed in [45| . The relevant part of 9\ is 



dV Trcr + / dS Tr (03(7757™ + C5CT75 + c 6 a) > -I , (E6) 



where the first integral is over the volume, while the second is over the boundary. The dots indicate terms that either 
vanish or do not depend on a and 7 m indicates the gamma matrices normal to the boundary. The coefficients C3, C5, 
C6 are universal multipliers and explicit forms for the case in question can be found in Ref. (45j . Restricting to the 
case of bag boundary conditions that are a special case of chiral boundary condition with the chiral angle set to zero, 
results of Ref. [45[ show that C5, cq > and C3 = 0. Thus, the overall sign of 9\ does not change, and the previous 
conclusions remain valid. Changing the boundary conditions allowing for a non-vanishing chiral angle modifies the 
coefficients c. However, for the case of fiat boundaries, the sign of 9\ remains positive. 

2. Large Condensates 

If the values of the condensate is large, one may proceed in a similar way to the high temperature case discussed 
above, owing to the presence of the exponential in (|E2[) that suppresses contributions from the large-i region of the 
integration range. Similar steps give 
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2(4^72^4 - r —° 1 / 2(J ~° ia - Z V ™ 3 /2^ + V ^ 5/2 - + (7 3 ^ + —^7/2^3- 

+0^ - log(ftr) 2 f^9 a 4 - a 2 9 Y + § 2 \\ } + °— . (E7) 



where for D = 3 



c = 46*o 8Li b/2 

ci = 2%/2# 1/2 SLi 2 
15 
2/32 
2^2, 
T 2 

105 * ._. 3 
32^° 4/3-2 

In the above expression we have defined the following combination of the polylogarithm functions 



c 2 = ^2~9o 5Li 7/2 + 29i 5Li 3/2 



2\> Z » /— - 

c.3 = -^-#1/2 SLi 3 + V 26*3/2 SLh 



c ± = ^r^°o SLi g/2 + — 9i SLi 5/2 + 9 2 SLi 1/2 



5Li a = Li a f- e -^+^ + Li a (-e-P^-^A (E8) 
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